SEPARATION  CONTROL  USING  ZNMF  DEVICES: 
FLOW  PHYSICS  AND  SCALING  LAWS 


FINAL  REPORT 

AFOSR  GRANT  FA9550-05- 1-0093 


By 

Ye  Tian  and  Louis  N.  Cattafesta  III 


Department  of  Mechanical  and  Aerospace  Engineering 
University  of  Florida 


20080604029 


1 


REPORT  DOCUMENTATION  PAGE  SKSSSmi.. 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  data  sources, 

gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection 
of  information,  including  suggestions  for  reducing  this  burden  to  Washington  Headquarters  Service,  Directorate  for  Information  Operations  and  Reports, 

1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget, 

Paperwork  Reduction  Project  (0704-0188)  Washington,  DC  20503. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1.  REPORT  DATE  (DD-MM-YYYY) 

2.  REPORT  TYPE 

Final  Technical  Report 

3.  DATES  COVERED  (From  -  To) 

15  Jan  2005-31  Dec  2007 

4.  TITLE  AND  SUBTITLE 

Separation  Control  Using  ZNMF  Devices:  Flow  Physics  and  Scaling 
Laws 

5a.  CONTRACT  NUMBER 

FA9550-05-1  -0093 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Dr.  Louis  N.  Cattafesta 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Department  of  Mechanical  and  Aerospace  Engineering 

University  of  Florida 

Gainesville  FL  32611 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Office  of  Scientific  Research  (AFOSR) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 
AFOSR 

875  N.  Arlington  St.,  Rm.  3112 

11.  SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

12.  DISTRIBUTION  AVAILABILITY  STATEMENT  .  r 

AF 

DISTRIBUTION  A:  Approved  for  public  release;  distribution  unlimited. 

rRL-SR-AR-TR-08-0293 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 


The  primary  goal  of  this  research  is  to  implement  a  closed-loop  control  system  to  control  separated  flow  and  to 
evaluate  the  performance  of  the  controller.  A  control  system  that  includes  an  array  of  actuators,  sensors  (pressure 
sensors  or  lift/drag  balance)  and  a  digital  controller  is  proposed  to  control  flow  separation  in  a  closed-loop  fashion.  The 
first  chapter  introduces  the  flow  physics  and  active  control  approaches  of  flow  separation.  It  is  organized  as  follows. 
First,  a  brief  overview  of  separation  control  is  provided  to  orient  the  reader,  followed  by  the  motivation.  Than  a  technical 
background  section  is  presented  to  review  previous  work  reported  in  the  literature.  Finally,  the  objectives  and  technical 
approaches  of  this  research  are  presented. 


15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF  RESPONSIBLE  PERSON 

ABSTRACT 

OF  PAGES 

a.  REPORT 

Unclassified 

b.  ABSTRACT 

Unclassified 

C.  THIS  PAGE 

Unclassified 

Unclassified 

104 

19b.  TELEPONE  NUMBER  ( Include  area  code) 

(703) 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI-Std  Z39-18 


Nomenclature 


.4 


Abbreviations . 5 

Abstract . 6 

1  Introduction . 7 

1 . 1  Overview . 7 

1.2  Motivation . 8 

1.3  Background . 8 

1 .3. 1  Two-Dimensional  Separation  Flow  Physics . 8 

1 .3.2  Effects  of  Flow  Separation . 9 

1 .3.3  Control  of  Flow  Separation . 9 

1 .3.4  Closed-Loop  Control  Algorithms . 1 5 

1.4  Objectives . 18 

1.5  Approach . 18 

1 .6  Outline  of  This  Dissertation . 1 8 

2  Theoretical  Background . 22 

2. 1  Optimization  Algorithms . 22 

2.1.1  Downhill  Simplex  Algorithm . 22 

2. 1 .2  Extremum  Seeking  Algorithm . 23 

2.2  System  Identification  Algorithms . 23 

2.2.1  ARMARKOV/LS  Algorithm . 25 

2.2.2  ARMARKOV/LS/ERA  Algorithm . 25 

2.2.3  Recursive  ARMARKOV/Tocplitz  Algorithm . 28 

2.3  Adaptive  Disturbance  Rejection  Algorithms . 28 

2.3.1  ARMARKOV  Disturbance  Rejection  Algorithm . 30 

3  Simulation  and  Validation  experiments . 36 

3. 1  Optimization  Simulations . 36 

3.1.1  Downhill  Simplex  Simulation  Results . 36 

3. 1 .2  Extremum  Seeking  Simulation  Results . 36 

3.2  Vibration  Control  Testbed  Setup . 36 

3.3  Results  of  the  Vibration  Control  Tests . 39 

3.3. 1  Computational  Tests . 39 

3.3.2  System  Identification . 39 

3.3.3  Adaptive  Disturbance  Rejection . 40 

4  experimental  setup  and  data  analysis  method . 62 

4. 1  NACA  0025  Airfoil  Model . 62 

4.2  Synthetic  Jet  Actuators . 62 


2 


4.3  Experimental  Methods . 62 

4.3.1  Flow  Visualization . 63 

4.3.2  Lift/Drag  Balance . 63 

4.3.3  Dynamic  Pressure  Transducers . 63 

4.3.4  Hot  Wire  Anemometry . 64 

4.4  Control  System  Hardware  and  Software . 64 

4.5  Higher  Order  Statistical  Analysis  (HOSA) . 65 

5  Results  and  Discussion . 7.3 

5.1  Dynamic  Feedback  Control . 73 

5.1.1  Experimental  Configuration . 73 

5. 1 .2  System  Identification . 73 

5.1.3  Disturbance  Rejection . 75 

5.2  Nonlinear  Control . 78 

5.2.1  Experimental  Configuration . 78 

5.2.2  Flow  Instabilities . 78 

5.2.3  Actuator  Calibration . 78 

5.2.4  Adaptive  Control  Results . 80 

6  Summary  and  Future  Work . 98 

7  List  of  References . I  (X) 


3 


CD 

CL 

cP 


CM 

D 

d(xjc ) 
fc 

fe 

fm 


wake 


F+ 

h 

J 


W 


L 

P 

P. 

q 


"j 

u. 

X 


SEP 


'  TE 


Nomenclature 

Airfoil  chord  length 
Drag  coefficient  ( D/qc ) 

Lift  coefficient  ( L/qc ) 

Static  pressure  coefficient  (— — — ) 

q 

Steady  momentum  coefficient  (J/qc) 

Oscillatory  momentum  coefficient  ( ( J)/qc ) 
Drag 

Pressure  recovery  coefficient 

Filter  cutoff  frequency 
Excitation  frequency 
Modulation  frequency 

Shedding  frequency  of  separated  flow  ( U„/XTE  ) 
Wake  shedding  frequency  (U„/Wwakf ) 

Reduced  excitation  frequency  ( fXTE/U„  ) 

Slot  width 

Steady  jet  momentum  ( pU^h ) 

Oscillatory  jet  momentum  ( pu^h ) 

Lift 

Static  local  pressure 

Free  stream  pressure 

Free  stream  dynamic  pressure  ( pU^/l) 

Oscillatory  jet  velocity 

Mean  jet  velocity 

Free  stream  velocity 

Distance  from  separation  points  to  trailing  edge 
Distance  from  excitation  slot  to  trailing  edge 


p  Air  density 

6  Boundary  layer  momentum  thickness 

o]  Variance  of  signal 

cr  Variance  of  noise 
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Abbreviations 


AM 

Amplitude  Modulation 

AOA 

Angle  Of  Attack 

BM 

Burst  Modulation 

DAQ 

Data  AcQuisition 

DSP 

Digital  Signal  Processing 

ID 

IDentification 

LDV 

Laser  Doppler  Velocimetry 

MSE 

Mean  Square  Error 

PIV 

Particle  Image  Velocimetry 

PM 
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Power  Spectral  Density 
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Abstract 

Flow  separation  has  severe  adverse  effects  on  performance  of  flow-related  devices  (e.g., 
lift  loss  of  aircrafts).  Active  control  of  separated  flow  has  received  extensive  attention  as  it  is 
able  to  mitigate  or  eliminate  flow  separation  effectively.  Most  research  has  been  open-loop  in 
nature  (i.e.,  manually  adjusting  control  inputs  to  achieve  best  results).  Closed-loop  control  of 
separated  flow  has  many  potential  advantages  over  open-loop  control,  namely  optimization  in 
multi-dimensional  domain  with  constraints,  adaptability  to  changing  flow  conditions,  etc.  In  this 
research,  adaptive  closed-loop  control  is  used  to  reattach  the  separated  flow  over  a  NACA  0025 
airfoil  using  multiple  zero-net-mass-flux  (ZNMF)  actuators  that  cover  the  central  33%  of  the 
airfoil  span.  In  particular,  two  distinct  approaches  are  used.  Adaptive  disturbance  rejection 
algorithms  are  used  to  apply  dynamic  feedback  control  of  separated  flow.  The  closed-loop 
control  results  show  ~  7  x  improvements  in  the  lift/drag  ratio,  with  a  corresponding  increase  in 
lift  and  reduced  drag  and  concomitant  reductions  in  the  fluctuating  surface  pressure  spectra.  On 
the  other  hand,  a  simplex  optimization  approach  uses  the  lift  and  drag  measured  by  a  strain- 
gauge  balance  for  feedback  and  searches  for  the  optimal  actuation  parameters  in  a  closed-loop 
fashion.  The  constrained  optimization  results  seeking  to  maximize  lift-to-drag  ratio  are 
promising  and  reveal  the  importance  of  forcing  nonlinear  interactions  between  the  shear  layer 
and  wake  instabilities. 
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1  Introduction 


The  primary  goal  of  this  research  is  to  implement  a  closed-loop  control  system  to  control 
separated  flow  and  to  evaluate  the  performance  of  the  controller.  A  control  system  that  includes 
an  array  of  actuators,  sensors  (pressure  sensors  or  lift/drag  balance)  and  a  digital  controller  is 
proposed  to  control  flow  separation  in  a  closed-loop  fashion. 

This  first  chapter  introduces  the  flow  physics  and  active  control  approaches  of  flow 
separation.  It  is  organized  as  follows.  First,  a  brief  overview  of  separation  control  is  provided  to 
orient  the  reader,  followed  by  the  motivation.  Then  a  technical  background  section  is  presented 
to  review  previous  work  reported  in  the  literature.  Finally,  the  objectives  and  technical 
approaches  of  this  research  are  presented. 

1.1  Overview 

Flow  separation  is  identified  as  one  of  the  most  important  flow  phenomena  due  to  its 
severe  adverse  effects  on  flow-related  devices.  Following  the  introduction  of  the  concept  of  the 
boundary  layer  by  Prandtl  (1904),  flow  separation  has  received  considerable  attention  in  the  fluid 
dynamics  community. 

Flow  separation  is  the  breakaway  or  detachment  of  fluid  from  a  solid  surface  (Greenblatt 
and  Wygnanski  2000).  Flow  separation  incurs  a  large  amount  of  energy/lift  loss  and  limits  the 
performance  of  many  flow-related  devices  (e.g.,  airplanes,  diffusers,  etc.).  Researchers  have 
been  trying  to  eliminate  or  at  least  mitigate  flow  separation  for  over  a  century  because  of  its  large 
potential  payoff  in  many  applications. 

As  shown  in  Figure  1-4,  control  of  separated  flow  is  divided  into  two  main  categories: 
active  control  and  passive  control.  Active  control  provides  external  energy  into  the  flow  while 
passive  control  does  not.  Some  passive  separation  control  methods,  such  as  geometrical  shaping 
and  turbulators  (i.e.,  turbulence  generators),  are  commonly  used  because  of  their  simplicity  and 
feasibility.  On  the  other  hand,  tremendous  progress  has  been  made  in  active  separation  control 
over  the  past  twenty  years.  Traditional  active  separation  control  methods,  such  as  steady 
blowing  and  suction,  were  initially  used  to  control  flow  separation  (Gad-el-Hak  2000).  These 
methods  were  able  to  control  of  separation  to  some  extent.  However,  they  were  far  from 
optimal  because  the  overall  energy  required  input  required  to  gain  a  meaningful  lift  increase  or 
drag  reduction  was  comparable  to  the  energy  saved  via  control  of  separation  (Greenblatt  and 
Wygnanski  2000). 

Schubauer  and  Skramstad  (1948)  first  introduced  a  breakthrough  in  active  flow  control: 
periodic  excitation.  This  technique  requires  much  less  energy  than  traditional  steady  active 
methods  and  accelerates  and  regulates  the  generation  of  large  coherent  structures  that  are 
primarily  responsible  for  the  transport  of  momentum  across  the  flow  (Greenblatt  and  Wygnanski 
2000).  The  increased  large  coherent  structures  make  the  flow  more  resistant  to  separation. 
Periodic  excitation  has  subsequently  been  shown  to  be  superior  to  steady  boundary  layer  control 
methods  by  many  researchers  (Seifert  1996;  Greenblatt  and  Wygnanski  2000;  Nishri  and 
Wygnanski  1998).  Because  of  these  reasons,  periodic  excitation  is  now  widely  used  to  control 
flow  separation.  Optimal  excitation  locations,  waveforms  shapes,  and  frequencies  of  periodic 
perturbations  have  been  systematically  studied  by  numerous  researchers  (Seifert  and  Pack 
2003 A,  Amitay  et  al.  2001).  Yet  none  of  these  studies  has  used  feedback  control  to  “optimize” 
the  excitation  waveform. 
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One  of  the  most  important  aspects  of  separation  control  is  the  actuation  mechanism  that 
introduces  periodic  perturbations  into  the  flow  structure.  Internal  acoustic  excitation  (Hsiao  et 
al.  1990;  Huang  et  al.  1987),  speakers  (Narayanan  and  Banaszuk  2003),  oscillatory  blowing 
valves  (Allen  et  al.  2000),  and  MEMS-based  actuators  (Rathnasingham  and  Breuer  2003),  etc. 
have  been  investigated.  Among  these,  synthetic  jet  or  zero-net  mass  flux  (ZNMF)  actuators  have 
been  the  focus  of  significant  research  for  the  past  decade  due  to  their  utility  in  flow  control 
applications  (Glezer  and  Amitay  2002).  ZNMF  actuators  utilize  the  working  fluid  and  do  not 
require  an  external  fluid  source,  which  makes  them  very  attractive  from  a  systems 
implementation  perspective.  Significant  progress  has  been  made  in  the  modeling  and  design  of 
such  devices  (Gallas  et  al.  2003,  2005).  More  details  of  the  synthetic  jet  actuators  used  in  this 
research  are  described  in  Chapter  4.  The  driving  frequency,  location,  and  momentum  coefficient 
of  the  actuation  are  the  primary  parameters  that  characterize  their  performance  (Amitay  et  al. 
2001). 

Although  separation  control  has  received  extensive  attention,  to  date  most  studies  have 
focused  on  open-loop  separation  control.  In  the  author’s  opinion,  this  open-loop  approach  is  due 
to  a  fluid  mechanics  bias  to  avoid  using  a  more  complex  closed-loop  control  approach.  Closed- 
loop  separation  control  has  the  potential  to  save  more  energy  than  open-loop  methods  (Cattafesta 
et  al.  1997)  and  make  separation  control  systems  adaptable  to  different  flow  conditions.  Few 
experimental  studies  have  focused  on  closed-loop  separation  control.  For  example,  Allan  et  al. 
(2000)  attempted  to  tune  a  PID  controller  for  closed-loop  separation  control  and  showed  that  the 
integral  gain  was  the  most  effective  as  a  result  of  the  large  time  constant  of  their  low  bandwidth 
actuator  system.  However,  the  realized  model  and  controller  were  simple.  Their  results  merely 
scratched  the  surface  of  what  can  possibly  be  accomplished.  Therefore,  it  is  believed  that  control 
of  flow  separation  using  an  array  of  high  bandwidth  actuators  and  surface  sensors  (pressure  or 
shear  stress)  is  an  excellent  candidate  for  closed-loop  separation  control.  Hence,  implementation 
of  feedback  controllers  including  more  advanced  modeling  and  control  algorithms  to  flow 
separation  control  is  proposed  and  is  the  focus  of  this  research. 

1.2  Motivation 

Numerous  applications  of  separation  control,  each  with  significant  potential  payoffs,  have 
been  identified  (Greenblatt  and  Wygnanski  2000).  Many  separation  control  strategies  have  been 
applied  on  civil  and  military  aircrafts  and  underwater  vehicles.  However,  most  of  the 
applications  are  open-loop  in  nature  because  of  their  simplicity.  Although  some  closed-loop 
separation  control  research  has  been  done  (Allan  et  al.  2000;  Banaszuk  et  al.  2003,  etc),  they  are 
not  sufficiently  developed  to  be  implemented  on  real  vehicles.  The  goal  of  this  research  is  to 
design  and  implement  various  closed-loop  control  systems  for  control  of  separated  flows  and  to 
seek  physical  insights  behind  the  control  schemes.  The  main  advantages  of  closed-loop 
separation  control  potentially  include  better  performance,  energy  efficiency  and  adaptability  to 
changing  of  flow  conditions. 

1.3  Background 

1.3.1  Two-Dimensional  Separation  Flow  Physics 

Under  the  circumstances  of  an  adverse  pressure  gradient  (dp/dx>0),  fluid  particles  are 
retarded  by  both  the  increasing  pressure  as  well  as  wall  skin  friction.  If  the  adverse  pressure 
gradient  is  of  sufficient  strength,  fluid  particles  near  the  wall  are  likely  to  separate  from  the  wall 
and  move  upstream.  This  is  due  to  the  fact  that  these  particles  have  finite  kinetic  energy  and 
cannot  penetrate  far  into  the  adverse  pressure  gradient  region.  The  flow  separates  from  the 
boundary  layer  and  forms  large  scale  vortical  structures  in  the  separated  region  (Figure  1-1). 
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Assuming  two-dimensional,  incompressible,  steady  flow  with  negligible  gravity,  the 
streamwise  (“  x  ”)  component  of  the  momentum  equation  at  the  wall  reduces  to 


or 

d2 u  1  dp 

v^=—r  (2) 

ay  p  dx 

where  v=p/p  is  the  kinematic  viscosity,  y  is  the  wall  normal  coordinate,  x  is  the  streamwise 
coordinate  with  a  corresponding  u  velocity  and  streamwise  pressure  gradient  dp/dx  . 

From  eqn.  (2),  we  can  see  that  only  an  adverse  pressure  gradient  ( dp/dx  >  0 )  can  cause  a 
point  of  inflection  in  the  velocity  profile  and  the  curvature  changing  sign  to  make  the  profile  S- 
shape.  In  this  case,  separation  will  occur  when  the  adverse  pressure  gradient  is  strong  enough  to 
make  the  right  hand  side  of  eqn.  (2)  positive  (shown  in  Figure  1-1). 

1.3.2  Effects  of  Flow  Separation 

In  the  separation  region,  the  normal  velocity  component  significantly  increases  as  well  as 
the  thickness  of  boundary  layer.  Therefore,  the  boundary  layer  approximations  are  no  longer 
valid  and  the  problem  can  no  longer  be  solved  using  boundary  layer  theory. 

Flow  separation  significantly  changes  the  pressure  distribution  around  the  surface.  Such 
deviations  are  usually  detrimental.  As  an  example.  Figure  1-3  shows  the  CL  and  CD  of  a 
NACA0025  airfoil  versus  angle  of  attack  measured  by  a  lift/drag  balance  at  Re  =  100,000. 
When  the  angle  of  attack  increases  from  zero  degree,  both  CL  and  CD  increase  as  expected. 
However,  CL  drops  dramatically  due  to  flow  separation  at  about  13  degrees  of  angle  of  attack. 
At  the  same  time,  CD  continues  to  increase  beyond  the  inception  of  stall.  Both  of  these  effects 

generally  have  a  negative  impact  on  the  airplane  performance.  However,  some  applications 
utilize  flow  separation.  For  example,  the  use  of  spoilers  on  airplanes  during  landing  reduces  the 
lift  and  increases  drag  to  allow  the  brakes  to  work  more  efficiently. 

More  commonly,  we  want  to  mitigate  or  eliminate  flow  separation.  Typical  applications 
of  flow  separation  control  include:  separation  control  of  various  airfoils  to  increase  CLmaj  for 

larger  payload  (Greenblatt  and  Wygnanski  2000;  Seifert  and  Pack  2002;  etc);  to  reduce  engine 
power  and  noise  at  takeoff  (Gad-el-Hak  2000);  to  increase  efficiency  of  diffusers  (i.e.  pressure 
recovery)  (Banaszuk  et  al.  2003);  etc. 

1.3.3  Control  of  Flow  Separation 

Because  of  the  effects  mentioned  above  and  the  large  potential  payoff,  researchers  have 
been  preoccupied  with  delaying  flow  separation  or  eliminating  it  entirely.  As  suggested  by 
Cattafesta  et  al.  (2003),  the  classification  of  flow  control  is  chosen  as  shown  in  Figure  1-4  to  be 
consistent  with  terminology  used  in  active  noise  and  vibration  control.  Active  control  is 
subdivided  into  open-loop  versus  closed-loop  control.  Closed-loop  control  can  be  further 
classified  into  quasi-static  versus  dynamic,  the  distinction  between  the  two  being  whether  or  not 
the  feedback  control  is  performed  on  a  time  scale  with  the  dynamical  scales  of  the  flow.  Since 
fluid  flows  are  inherently  nonlinear  (Wu  et  al.  1998),  the  standard  frequency  preservation  of  a 
linear  system  does  not  hold.  Consequently,  nonlinear  feedback  control  on  a  very  slow  time 
compared  to  the  characteristic  times  scales  of  the  flow  is,  in  fact,  possible  and  attractive.  In 
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essence,  this  so-called  quasi-static  control  becomes  a  nonlinear  optimization  problem.  This 
research  will  investigate  both  classes  of  closed-loop  control  shown  in  Figure  1-4. 

Other  fluid  dynamic  issues  have  been  studied  extensively,  such  as  the  effects  of  Reynolds 
number,  frequency,  actuator  and  sensor  locations,  momentum  coefficient,  surface  curvature,  and 
compressibility,  etc..  Although  the  topic  of  this  research  is  closed-loop  separation  control,  the 
results  and  conclusions  from  the  open-loop  control  studies  should  serve  as  a  sound  physical  basis 
for  effective  control  and  are  reviewed  below. 

1.3.3.1  Open-loop  separation  control. 

Periodic  excitation  has  been  shown  to  be  much  more  effective  than  steady  forcing 
because  it  enhances  the  momentum  transport  across  the  flow  domain  at  a  substantial  reduction  in 
energy  expenditure.  It  accelerates  and  regulates  the  generation  of  large  coherent  structures  that 
are  primarily  responsible  for  the  momentum  transport  across  the  flow  (Greenblatt  and 
Wygnanskj  2000).  The  enhanced  momentum  transport  forces  the  separated  flow  to  reattach  to 
the  surface  and  form  a  thick  turbulent  boundary  layer  in  a  time-averaged  sense.  The 
reattachment  of  the  boundary  layer  regains  the  pressure  suction  zone  on  the  upper  surface  of  the 
airfoil  and  thus  enhances  the  lift  performance.  Furthermore,  the  superposition  of  weak  suction 
on  the  periodic  excitation  enhances  the  receptivity  of  the  separated  shear  layer  to  the 
fundamental  excitation  frequency  and  thus  the  effectiveness  of  periodic  excitation  (Seifert  and 
Pack  2002). 

Given  the  improved  performance  of  periodic  excitation  to  control  flow  separation, 
researchers  have  sought  to  optimize  separation  control  via  time-consuming  parametric  variations. 
Significant  parameters  or  conditions  that  affect  the  performance  of  separation  control  have  been 
identified.  Although  they  are  discussed  separately  below,  one  should  keep  in  mind  that  these 
factors  are  all  coupled  with  each  other. 

Actuation  frequency.  First,  consider  the  characteristic  flow  structures  associated  with 
separated  flow.  Based  on  previous  studies,  Mittal  et  al.  (2005)  summarize  the  three  situations 
with  regards  to  separated  flow,  as  shown  in  Figure  1-5.  In  post-stall  flow  (case  C  in  Figure  1-5), 
leading-edge  shear  layer  rollup  and  vortex  shedding  in  the  wake  are  two  characteristic  features 
(Wu  et  al.  1998).  Huerre  and  Monkewitz  (1990)  suggest  that  this  type  of  shear  flow  (with  a 
pocket  of  absolute  instability  of  sufficient  size)  may  display  intrinsic  dynamics  of  the  same 
nature  as  in  a  closed-flow  system,  in  which  disturbances  can  grow  upstream  (i.e.  global 
instability).  Therefore,  it  is  reasonable  to  postulate  that  separated  flow  over  an  airfoil  acts  as  a 
nonlinear  multi-frequency  closed-flow  system.  In  such  a  system,  the  shear  layer  instability  (with 
characteristic  frequency  fSL)  and  the  global  wake  instability  (with  vortex  shedding  frequency 

fwake )  may  interact  with  each  other  in  a  nonlinear  fashion.  In  case  B,  a  closed  separation  bubble 
is  present  at  some  distance  downstream  of  the  leading  edge.  In  this  case  there  are  potentially 
three  characteristic  flow  frequencies  in  the  separated  flow:  fSL,  fwakeand  f  ,  where  the  new 

scale,  fscp ,  corresponds  to  the  characteristic  frequency  of  the  separation  bubble. 

The  scales  of  the  three  frequencies  are  fSL  ~  U/0SL  ,  fSL~U/Lsep  and  fwake~  U/Wwake , 
where  0SL  is  the  shear  layer  thickness,  Lscp  is  the  length  of  the  separation  bubble  and  Wwake  is 
the  width  of  wake.  Prasad  and  Williamson  (1996)  also  show  that  fSL=ARe8fwake ,  where 
A  =  0.0235  and  B  =  0.67  .  Since  there  are  different  relevant  length  scales  that  are  included  in 
the  three  characteristic  frequencies,  one  should  expect  a  significant  variation  in  the  observed 
frequency  scales  and  the  corresponding  optimal  frequency. 
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The  present  study  is  focused  on  how  flow  systems  respond  to  modulated  (e.g.  AM,  BM, 
PM)  unsteady  excitations  by  ZNMF  devices  targeting  the  inherent  flow  instabilities  that  lead  to 
the  presence  of  these  characteristic  flow  frequencies.  The  goal  is  to  search  for  optimal  forcing 
schemes  that  most  effectively  mitigate  flow  separation  via  nonlinear  interaction  of  the 
instabilities. 

Much  research  has  been  conducted  to  determine  what  excitation  frequencies  are  most 
effective  for  separation  control.  However,  except  for  the  general  agreement  that  periodic 
excitation  is  far  more  effective  than  steady  blowing,  the  range  of  optimal  actuation  frequencies  is 
a  current  subject  of  intense  debate.  A  dimensionless  actuation  frequency  is  typically  defined  for 
this  purpose.  However,  three  slightly  different  definitions  have  been  given  for  a  so-called 
dimensionless  frequency  F+ :  1)  F+  =  feXTE/Uoo ,  where  fe  is  the  excitation  frequency,  XTC  is 
the  distance  from  the  excitation  slot  to  the  trailing  edge  and  U„  is  the  free  stream  velocity;  2) 
F+  =  feLsep/Uoo ,  where  Lsep  is  the  distance  from  separation  to  reattachment;  and  3)  F+  =  fec/Uoo  , 

where  c  is  the  chord  length.  These  three  are  nearly  identical  for  post-stall  flow  (where  the 
separation  bubble  length  is  approximately  the  airfoil  chord),  but  they  scale  very  differently  if  a 
closed  separation  bubble  of  finite  extent  is  present.  One  should  notice  that  none  of  these 
definitions  is  related  to  the  shear  layer  frequency  (fSL).  Most  researchers  implicitly  ignore  this 

important  frequency  when  studying  separation  control. 

Herein,  some  results  regarding  actuation  frequency  in  previous  studies  are  summarized. 
Among  studies  that  define  F+  =  feXTE/U„  ,  Wygnanski  and  his  colleagues  conclude  that  the 

optimal  excitation  frequency  is  of  order  unity  F+=0(1)  (Seifert  et  al.  1996,  Nishri  and 
Wygnanski  1998,  Greenblatt  and  Wygnanski  2000)  and  have  found  that  so-called  high  frequency 
fording  F+=O(10)  is  ineffective  for  their  airfoil  (NACA  0015)  and  flow  conditions.  Conversely, 
using  the  same  definition  of  F+ ,  Amitay  et  al.  (2001)  found  that  when  the  excitation  frequency 
F+  >0(10),  the  lift-to-pressure  drag  ratio  was  larger  than  that  when  the  excitation  frequency 

F+  <  4  .  Honohan  et  al  (2000)  also  suggested  that  higher  reduced  frequencies  (F+  >10)  can  be 
effective.  They  argued  that  it  is  because  the  high  frequency  excitation  produces  a  virtual 
aerodynamic  surface  modification  that  thins  the  turbulent  boundary  layer  and  results  in  a  local 
favorable  pressure  gradient. 

Besides  this  argument,  there  may  be  two  other  possible  reasons  accounting  for  this 
interesting  discrepancy.  First,  the  length-scale  XTC  may  not  be  appropriate  for  their  airfoil 

because  of  the  formation  of  a  closed  separation  bubble.  Instead,  if  were  used,  this 
discrepancy  might  not  exist.  Second,  as  mentioned  earlier,  the  shear  layer  frequency  fSL  may 
also  be  important  (Mittal  et  al  2005).  Here,  fSL«U/0,  where  0  is  the  boundary  layer 
momentum  thickness  and  not  XTE  or  Lsep.  The  different  frequency  scales  are  indicative  of 

different  flow  instabilities  that  may  exist  in  the  flow  and,  if  present,  may  compete  with  each 
other  (Wu  et  al  1998).  When  periodic  excitation  is  introduced,  one  or  more  of  these  instabilities 
may  be  energized.  The  controlled  flow  may  then  be  regulated,  and  thus  lift  performance  may  be 
enhanced.  This  may  explain  the  observed  variations  of  the  optimal  excitation  frequency. 

Along  these  lines,  an  innovative  forcing  approach  that  uses  multiple  harmonically  related 
frequencies  is  presented  by  Narayanan  and  Banaszuk  (2003).  They  demonstrated  improvements 
of  this  new  approach  versus  single  frequency  sinusoidal  forcing  in  control  of  separation  in  a 
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diffuser,  although  its  effectiveness  requires  further  investigation.  To  extend  this  idea  further,  one 
can  use  excitations  with  multiple  frequency  components  corresponding  to  the  characteristic 
frequencies  mentioned  above.  This  idea  will  be  investigated  in  this  research. 

Excitation  amplitude.  Another  key  control  parameter  in  a  ZMNF  device  is  jet  velocity 
Vj  (some  characteristic  velocity  measure,  e.g.  the  peak  or  an  average  velocity).  In  the  literature, 

the  jet  frequency  is  usually  non-dimensionalized  as  F+=fLscp/Uoo,  where  Lsep  is,  for  example,  the 
length  of  separation  region  and  U„  is  the  free  stream  velocity.  The  jet  velocity  is  usually  non- 
dimensionalized  by  U,,,.  Various  researchers  have  shown  that  control  authority  varies 
monotonically  with  Vj/U„,  for  a  sinusoidal  excitation  up  to  some  maximum  value  (Seifert  et  al. 

1993,  1996,  1999;  Glezer  and  Amitay  2002;  Mittal  and  Rampunggoon  2002).  In  practice, 
especially  in  high  speed  flows,  control  authority  is  often  lacking.  From  an  efficiency  standpoint, 
it  is  desirable  to  control  a  flow  with  minimal  actuator  input. 

Modulation  signals.  Piezoelectric  actuators  have  fast  dynamic  response  and  low  power 
consumption.  However,  the  use  of  piezoelectric  actuators  has  been  limited  because  of  the 
diminution  in  their  response  outside  a  narrow  frequency  band  around  their  resonance  frequency 
and  the  need  for  testing  over  a  wide  frequency  range  due  to  the  issues  discussed  in  the  last 
section. 

Wiltse  and  Glezer  (1993)  introduced  a  clever  amplitude  modulation  method  to  flow 
control  problems  to  overcome  this  problem.  The  piezoelectric  actuator  is  resonantly  driven  with 
a  carrier  waveform,  e(t) ,  which  is  amplitude  modulated  with  a  time-harmonic  wave  train: 

e(t)=[  1  +esin(o)m  t+cpm )]  Arsin(wct)  (3) 

where  Ar  is  the  amplitude  of  the  carrier  signal,  e  is  the  degree  of  modulation  (0<  e  <  1 ),  o>c  is 
the  carrier  frequency  (or  the  resonant  frequency  of  the  actuator)  in  rad/s ,  a)m  is  the  modulation 
frequency  (which  is  also  the  desired  excitation  frequency  or  receptive  frequency  of  the  flow)  in 
rad/s,  and  <pm  is  the  phase  of  the  modulating  signal.  By  using  trigonometric  identities,  one  can 

show  that  e(t)  contains  frequency  components  at  coc  and  o)c±o)m,  However,  when  the 
excitation  amplitude  is  high  enough,  e(t)  is  demodulated  by  the  nonlinear  fluid  dynamical 
system  that  is  associated  with  the  formation  and  coalescence  of  nominally  spanwise  vortices. 
This  nonlinearity  results  in  the  presence  of  coc  and  toc±tom  and  also  (0m  in  the  flow.  In  practice, 
coc  is  set  at  the  resonance  frequency  of  the  piezoelectric  actuator  (which  is  usually  »com)  and 
com  is  set  at  the  desired  low  frequency  corresponding  to  the  desired  excitation  frequency  fe . 

Along  these  lines,  other  modulation  signals  such  as  burst  modulation  and  pulse 
modulation  can  also  be  used.  This  modulation  technique  allows  the  actuator  operating  at  its 
resonant  frequency  to  generate  a  significant  flow  disturbance  while  effectively  manipulating 
flows  at  characteristic  frequencies  of  the  flow.  It  provides  a  much  more  flexible  approach  than 
matching  the  resonant  frequency  of  the  actuator  with  the  receptive  frequencies  of  the  flow. 

However,  some  features  of  the  technique  should  also  be  kept  in  mind.  First,  the  actuator 
is  driven  continuously  near  its  resonant  frequency,  so  the  probability  of  mechanical  failure  is 
greater  than  when  it  is  driven  off  resonance.  Second,  as  mentioned  above,  demodulation  of  the 
waveform  is  due  to  nonlinearities  of  the  flow  and  actuator.  As  a  result,  feedback  controllers 
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designed  based  on  a  linear  assumption  may  not  work  as  desired.  This  aspect  will  be  studies  in 
this  research. 

Actuation  location.  It  is  argued  by  many  researchers  that  the  optimal  actuation  location 
is  at  the  vicinity  of  the  point  of  separation  (Amitay  et  al.  2001,  Seifert  et  al.  1996,  Seifert  and 
Pack  2003).  This  is  physically  plausible  since  the  disturbances  introduced  at  this  location  can 
most  effectively  transport  momentum  between  the  free  shear  layer  and  the  separated  region. 
However,  this  has  not  been  systematically  studied  because  of  some  practical  limitations,  namely 
the  difficulty  of  installing  multiple  actuators  inside  an  airfoil.  Amitay  et  al.  (2001)  used  an 
unconventional  airfoil  that  had  an  aft  portion  of  a  symmetric  airfoil  attached  to  a  circular 
cylinder  forebody  with  a  synthetic  jet  slot  that  could  be  adjusted  by  rotating  the  cylinder.  They 
state  that  the  closer  the  control  is  located  to  the  observed  separation  point,  the  less  power  is 
required  to  reattach  the  flow.  They  also  made  an  interesting  point  that  if  either  the  separation 
location  is  unknown  or  practical  limitations  preclude  control  near  the  separation  location,  the 
momentum  coefficient  CM  may  be  manipulated  to  achieve  optimal  performance. 

Besides  the  effects  of  actuation  location  discussed  above,  the  interaction  of  adjacent 
synthetic  jet  actuators  has  been  investigated  by  Holman  et  al.  (2003).  They  found  that  relative 
phasing  between  adjacent  actuators  does  not  appear  to  affect  the  effectiveness  of  separation 
control  significantly  for  their  airfoil  (NACA  0025)  and  flow  conditions  (R,=105  and 

AO  A  =  12°). 

In  summary,  based  on  the  previous  studies  it  is  suggested  that  slightly  upstream  of  the 
separation  location  is  the  “best”  place  to  introduce  actuation.  Furthermore,  a  combination  of 
upstream  leading  edge  and  downstream  trailing  edge  actuations  may  also  be  a  good  candidate 
and  remains  to  be  investigated  (Mittal  et  al  2005).  Wu  et  al  (1998)  discuss  this  idea  in  the 
context  of  the  Kutta-Joukowski  lift  formula  (L=-pUT),  which  assumes  the  flow  is 
incompressible  and  steady.  In  the  formula,  L  is  the  lift,  U  is  free  stream  velocity  and  T  is  the 
circulation  (a  counterclockwise  circulation  is  assumed  positive).  Although  the  separation  is  an 
unsteady  process,  this  formula  still  holds  in  a  time-averaged  sense  for  the  entire  flow.  Based  on 
these  arguments,  if  the  combination  of  leading  edge  and  trailing  edge  actuation  can  be  designed 
to  alter  the  circulation  of  the  airfoil,  it  should  be  able  to  control  flow  separation  in  some  manner. 

Effects  of  Reynolds  number  and  compressibility.  It  is  shown  that  control  of  flow 
separation  is  insensitive  to  the  Reynolds  number  at  high  chord  Reynolds  numbers  of  11-30 
million  (Seifert  and  Pack  2003  A,  B,  Greenblatt  and  Wygnanski  2000).  The  Reynolds  number 
has  a  very  weak  effect  on  pressure  distributions  around  the  surface,  regardless  of  the  Mach 
number. 

On  the  other  hand,  strong  Reynolds  number  effects  are  identified  in  the  airfoil  baseline 
performance  at  moderately  compressible  flow  conditions  (Seifert  and  Pack  2001).  Reynolds 
number  effects  weaken  as  the  Mach  number  increases  and  a  stronger  shock  wave  develops. 
Compressibility  tends  to  elongate  the  separation  bubble  and  reduce  the  capability  of  periodic 
excitation  to  shorten  the  separation  bubble  with  similar  excitation  frequencies  and  momentum 
(Seifert  and  Pack  2001). 

It  is  also  suggested  by  Seifert  and  Pack  (2001)  that  in  the  presence  of  shock  waves  the 
excitation  location  should  be  slightly  upstream  of  the  shock  wave.  If  the  excitation  is  introduced 
well  upstream  of  the  shock  wave,  it  has  a  detrimental  effect  on  lift,  drag  and  wake  steadiness. 
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1. 3.3.2  Closed-loop  separation  control 

Closed-loop  experimental  separation  control  has  not  yet  received  significant  attention. 
This  section  first  reviews  some  development  of  the  micro-electro-mechanical  systems  (MEMS) 
based  actuators  because  of  their  potential  importance  to  high  bandwidth  closed-loop  control 
systems.  Then  the  limited  previous  work  on  closed-loop  separation  control  is  presented. 

For  closed-loop  flow  control  systems,  the  desired  actuators  should  be  fast,  power 
efficient,  and  reliable.  In  previous  separation  control  studies,  acoustic  excitation  (Hsiao  et  al. 
1990  and  Huang  et  al.  1987)  seems  facility  dependent  because  the  acoustic  drivers  stimulate  the 
wind  tunnel  resonant  modes  to  excite  the  separated  flow;  oscillatory  blowing  valves  (Allen  et  al. 
2000)  appear  to  have  slow  dynamic  response;  active  flexible  wall  transducers  (Sinha  2001)  have 
complicated  structures  despite  its  high  actuation  efficiency  and  ability  to  actuate  and  sense  with 
the  same  hardware.  These  drawbacks  have  limited  the  use  of  these  actuators. 

On  the  other  hand,  synthetic  jet  (ZNMF)  actuators  have  been  the  focus  of  significant 
research  activities  for  the  past  decade  due  to  their  utility  in  flow  control  applications  (Glezer  and 
Amitay  2002).  They  utilize  the  working  fluid  and  do  not  need  external  fluid  injection.  They  can 
force  the  momentum  transfer  across  the  flow  without  net  mass  flux  (thus  the  name  “synthetic”). 
The  design  of  synthetic  jets  is  also  flexible  and  the  working  frequency  range  can  be  tuned 
according  to  different  flow  control  applications.  In  addition,  the  recent  paper  by  Gallas  et  al. 
(2003)  presents  a  lumped  element  model  of  a  piezoelectric-driven  synthetic  jet  actuator.  They 
provide  a  novel  method  to  design  and  model  synthetic  jets,  which  makes  them  very  suitable  for 
closed-loop  separation  control.  In  lumped  element  modeling  (LEM),  the  individual  components 
of  a  synthetic  jet  are  modeled  as  elements  of  an  equivalent  electrical  circuit  using  conjugate 
power  variables  (i.e.,  power  =  generalized  flow  x  generalized  effort).  The  frequency  response 
function  of  the  circuit  is  derived  to  obtain  an  expression  for  Qoul/VAC  ,  the  volume  flow  rate  per 

applied  voltage.  The  comparison  between  the  LEM  and  experimental  frequency  response  is 
shown  in  Figure  1-6. 

For  a  variety  of  reasons,  closed-loop  control  in  a  real-time  experiment  has  been 
traditionally  difficult  to  achieve.  In  reduced-scale  laboratory  experiments,  the  characteristic 
frequencies  of  separated  turbulent  flows  are  proportionally  higher  than  those  on  full-scale 
models,  which  requires  high  frequency  sensing  and  actuating  capabilities.  Furthermore,  real¬ 
time  experiments  require  the  digital  control  system  to  sample  at  a  minimum  of  twice  of  the 
highest  frequency  of  interest.  The  availability  of  hardware  (including  actuators,  sensors  and  real¬ 
time  control  systems)  therefore  imposes  significant  limitations  on  the  complexity  of  the  closed- 
loop  control  system.  Lower  order  system  models  are  typically  required  to  reduce  the  complexity 
of  the  system. 

Many  model-based  approaches  are  being  developed  and  have  shown  promising  results. 
Proper  Orthogonal  Decomposition  (POD)  based  low  order  models  have  been  studied  extensively 
(Holmes  et  al.  1998;  Tadmor  et  al.  2007)  owing  to  their  relatively  high  resolution  and  low 
computational  intensity.  Other  reduced-basis  models  have  also  been  studied  (Coller  et  al.  2000; 
Wang  et  al.  2003).  These  models  require  that  multiple  measurements  are  simultaneously 
available  in  the  flow  field.  However,  this  is  impractical  in  feedback  separation  control  and 
surface  measurements  are  required  in  most  applications.  Ausseur  et  al.  (2007)  implemented  a 
POD/mLSM  proportional  feedback  control  using  the  velocity  field  and  surface  pressure  data  to 
delay  flow  separation. 

Some  non-model  based  control  approaches  have  gained  favor  because  they  bypass  the 
complication  of  modeling  separated  flow  while  focusing  on  the  primary  control  objectives.  For 
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example,  Banaszuk  et  al.  (2003)  and  Becker  et  al.  (2006)  used  an  extremum-seeking  closed-loop 
control  algorithm  to  optimize  the  pressure  recovery  and  lift,  respectively.  The  present  author  in 
Tian  et  al.  (2006)  used  a  multi-dimensional  optimization  algorithm  to  optimize  lift-to-drag  ratio 
over  an  airfoil.  These  approaches  are  capable  of  “training”  the  excitation  signals  to  be  most 
effective  in  terms  of  the  objective  functions  (i.e.,  pressure  recovery,  lift-to-drag  ratio,  etc.).  The 
main  drawback  of  the  above  approaches  is  that  they  operate  on  a  time  scale  that  is  much  larger 
than  that  of  the  flow  dynamics.  In  other  words,  they  work  on  time-averaged  objective  functions 
by  explicitly  taking  advantage  of  the  nonlinear  nature  of  the  fluid  dynamics.  This  approach  has 
the  drawback  of  having  to  deal  or  cope  with  the  nonlinear  dynamics  with  no  guarantee  of 
success.  This  kind  of  approach  is  an  example  of  the  quasi-static  control  scheme  shown  in  Figure 
1-4. 

On  the  other  hand,  the  dynamic  feedback  control  is  used  to  model  and  control  separated 
flow  structures  based  on  surface  pressure  data  alone.  The  well-developed  adaptive  system 
identification  (ID)  algorithms  in  the  controls  community  are  utilized  to  model  the  flow  system 
dynamics  between  the  actuators  and  unsteady  surface  pressure  sensors.  The  system  ID 
algorithms  generate  known  actuation  signals  and  relate  these  signals  with  the  surface  pressure 
response  measured  by  sensors.  Linear  dynamical  equations  are  then  used  to  model  the 
relationship  in  a  gradient  descent  sense  (Haykin  2002).  The  system  therein  includes  the  dynamics 
of  the  actuators,  the  flow  structures  excited  by  the  actuation,  and  the  dynamics  of  the  sensors. 
The  system  information  is  then  used  to  predict  the  subsequent  evolution  of  the  pressure 
fluctuations.  Control  is  applied  using  a  spanwise  zero-net-mass-flux  (ZNMF)  actuator  slot  by 
attempting  to  reduce  the  power  of  the  surface  pressure  fluctuations  in  a  closed-loop  fashion,  thus 
suppressing  the  unsteady  flow  fluctuations  based  on  predicted  flow  characteristics.  A  similar 
idea  has  been  applied  to  control  of  flow-induced  cavity  oscillations  (Cattafesta  et  al.  1999)  and 
turbulent  boundary  layer  control  (Rathnasingham  and  Breuer  2003).  This  kind  of  approach  can 
be  categorized  as  a  dynamic  control  scheme  shown  in  Figure  1-4. 

1.3.4  Closed-Loop  Control  Algorithms 

According  to  the  classification  in  Figure  1-4,  the  control  algorithms  can  be  divided  into 
two  categories:  quasi-static  and  dynamic.  Optimization  algorithms  are  used  in  this  research  as 
quasi-static  algorithms.  They  are  used  to  optimize  target  functions  (such  as  lift,  pressure 
recovery,  etc.)  in  a  recursive  but  static  or  time-averaged  fashion.  On  the  other  hand,  recursive 
system  identification  and  disturbance  rejection  algorithms  are  widely  used  in  active  noise  control 
area  as  dynamic  algorithms.  No  one  has  attempted  to  apply  these  algorithms  to  the  closed-loop 
separation  control  problem.  This  section  gives  a  brief  review  of  the  two  types  of  the  algorithms. 
Details  will  be  given  in  chapter  2. 

1. 3.4.1  Optimization  algorithms 

Optimization  algorithms  are  widely  used  by  decision-makers  (e.g.  economists, 
governments).  They  often  need  to  choose  an  action  to  optimize  target  or  cost  functions,  such  as 
income,  profit,  etc.  In  a  typically  optimization  problem,  one  is  given  a  single  function  f  that 
depends  on  one  or  more  independent  variables.  The  goal  is  to  find  the  value  of  those  variables 
where  f  is  a  maximum  or  a  minimum  value.  In  this  research,  various  optimization  algorithms 
are  used  to  maximize/minimize  different  cost  functions,  such  as  lift,  drag  and  pressure  recovery. 
When  using  the  optimization  algorithms,  some  constraints  are  typically  included  in  the 
algorithms.  For  example,  one  often  seeks  to  limit  the  energy  expenditure  while  optimizing  the 
cost  function.  One  should  also  notice  that,  unlike  the  applications  used  by  the  decision-makers, 
the  cost  functions  used  in  this  research  are  measured  by  sensors  instead  of  analytical  functions. 
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Some  established  minimization  and  maximization  algorithms  are  summarized  by  Press  et 
al.  (1992).  Most  optimization  algorithms  can  be  easily  implemented  in  a  multi-dimensional 
space.  The  downhill  simplex  algorithm  and  the  Powell’s  algorithm  do  not  require  derivative 
calculations.  Between  these  two  algorithms,  the  downhill  simplex  algorithm  is  more  concise  and 
self-contained.  Both  of  them  require  storage  of  order  N2 ,  where  N  is  the  number  of  dimensions 
or  independent  variables.  Two  other  algorithms,  the  conjugate  gradient  and  quasi-Newton 
methods,  do  require  the  calculation  of  derivatives.  The  conjugate  gradient  method  requires  only 
order  N  storage,  while  the  quasi-Newton  method  requires  storage  of  order  N2.  On  the  other 
hand,  none  of  the  algorithms  mentioned  above  are  guaranteed  to  find  a  global  extremum.  They 
can  lead  to  local  extrema.  Finding  a  global  extremum  is  actually  a  very  difficult  problem.  Two 
standard  methods  are  typically  used  to  improve  the  probability  of  finding  a  global  extremum:  1) 
search  for  local  extrema  from  various  initial  conditions  and  pick  the  most  extreme  of  these;  2) 
perturb  a  local  extremum  to  see  if  the  algorithm  goes  back  to  the  same  value  or  finds  a  better 
result. 

There  are  several  global  search  algorithms  that  are  currently  active  in  research(e.g. 
Genetic  Algorithms  (GA)  (Holland  1975),  Particle  Swarm  Optimization  (PSO)  (Kennedy  1997) 
and  Simulated  Annealing  Method  (Haftka  and  Giirdal  1992)).  The  genetic  algorithms  and  the 
particle  swarm  optimization  are  both  derived  from  biology.  They  are  population-based 
algorithms,  namely  they  generate  a  population  of  points  at  each  iteration  and  the  population 
approaches  an  optimal  solution.  The  GA  and  PSO  take  advantage  of  the  large  search  population 
to  increase  probability  of  approaching  a  global  optimum.  The  simulated  annealing  method  is  an 
analogy  with  thermodynamics,  especially  with  the  way  metals  cool  and  anneal,  in  which  process 
nature  finds  the  minimum  energy  state.  The  essence  of  the  algorithm  is  to  allow  increase  of  cost 
function  with  some  probability  to  improve  the  changes  to  find  a  global  minimum. 

Another  optimization  algorithm  that  has  been  applied  to  flow  control  problems  is  called 
the  extremum-seeking  algorithm.  As  a  self-optimizing  control  algorithm,  the  extremum-seeking 
control  was  first  introduced  in  the  1950s.  After  Krstic  and  Wang  (1999)  provided  the  stability 
studies,  there  has  been  a  resurgence  of  interest  of  this  control  algorithm.  Banaszuk  et  al.  (2003) 
attempted  to  use  this  algorithm  in  the  diffuser  separation  control  problem.  They  were  successful 
in  maximizing  the  pressure  recovery  in  the  diffuser.  They  also  used  this  algorithm  to  control 
combustion  instability  (Banaszuk  et  al.  2000). 

1. 3.4.2  System  identification  and  disturbance  rejection  algorithms 

System  identification  and  disturbance  rejection  technologies  are  well  developed  and 
various  algorithms  are  available  in  the  active  noise  control  area.  Cattafesta  et  al.  (1999)  have 
applied  these  algorithms  to  other  flow  control  problems,  such  as  cavity  resonance  control.  No 
one  has  attempted  to  apply  this  kind  of  approach  to  the  separation  control  problem.  In  this 
research,  this  approach  is  investigated.  Some  system  identification  and  disturbance  rejection 
algorithms  are  reviewed  in  this  section. 

System  identification  algorithms.  In  general,  system  identification  (ID)  uses  measured 
signals  (i.e.,  inputs  and  outputs  of  the  system)  to  identify  (or  estimate)  the  unknown  system 
dynamics.  It  provides  necessary  system  information  for  control  algorithms.  System 
identification  algorithms  can  be  divided  into  two  categories:  offline  (or  batch)  and  online  (or 
recursive).  Offline  algorithms  first  acquire  data  and  then  try  to  estimate  a  low-order  dynamical 
system  model  using  these  data  offline.  Online  algorithms  identify  systems  recursively  while 
acquiring  data  in  real-time.  Online  system  identification  is  also  known  as  adaptive  filtering. 
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Least  square  (LS)  identification  algorithm  is  a  generally  used  offline  algorithm.  Akers 
and  Bernstein  (1997  A)  applied  this  approach  to  the  ARMARKOV/LS  identification  algorithm 
with  an  ARMARKOV  representation  (see  Chapter  2  for  a  detailed  description  of  the  algorithm). 
The  ARMARKOV/LS  identification  algorithm  uses  vectors  comprised  of  input-output  data  with 
a  least-squares  criterion  to  estimate  a  weight  matrix  containing  a  specified  number  of  Markov 
(i.e.,  pulse  response)  parameters  of  the  system.  Then  the  eigensystem  realization  algorithm 
(ERA)  (Juang  1994)  is  used  to  construct  a  minimal  state  space  realization  of  the  system.  This  is 
referred  to  as  the  ARMARKOV/LS/ERA  identification  algorithm. 

The  ARMARKOV/LS/ERA  identification  algorithm  has  two  clear  advantages  compared 
to  the  ARMA/LS  identification  algorithm  (Akers  and  Bernstein  1997  A).  First,  eigenvalues  of 
the  ARMARKOV  representation  are  less  sensitive  to  noise  compared  with  eigenvalues  of  the 
ARMA  representation.  Second,  the  singular  value  decomposition  of  a  block  Hankel  matrix 
constructed  from  the  estimated  Markov  parameters  provides  an  efficient  model  order  indicator 
(Juang  1994,  pp.  139). 

As  far  as  online  algorithms  are  concerned,  the  least-mean-square  (LMS)  algorithm  is  the 
most  commonly  used  algorithm.  A  more  computationally  intensive  algorithm  called  the 
recursive-least-square  (RLS)  algorithm  has  faster  convergence  and  smaller  steady-state  error 
than  the  LMS  algorithm  (Haykin  2002)  but  is  more  computationally  intensive.  Two  different 
types  of  structure  that  can  be  applied  to  each  of  the  algorithms  are  the  finite-impulse-response 
(FIR)  and  the  infinite-impulse-response  (HR)  filters.  The  FIR  filter  is  widely  used  due  to  its 
simple  architecture  and  inherent  stability  as  an  all-zero  model.  However,  its  simple  structure 
introduces  difficulties  for  a  system  with  low  damping.  The  IIR  filter  can  solve  this  problem  with 
significantly  lower-order  and,  therefore,  lead  to  reduced  computational  expense.  Unfortunately, 
the  disadvantages  of  an  IIR  filter  include  more  complicated  adaptive  algorithms  compared  with 
an  FIR  filter  and  the  possible  stability  problems  introduced  by  the  pole(s)  in  the  model  (Haykin 
2002;  Shynk  1989;  Netto  and  Diniz  1995). 

Applying  the  LMS  algorithm  to  the  ARMARKOV  representation,  Akers  and  Bernstein 
(1997  B)  introduced  the  recursive  ARMARKOV/Toeplitz  algorithm  that  is  based  upon  recursive 
identification  of  the  Markov  parameters  of  a  system.  It  estimates  the  Markov  parameters 
recursively  using  time-domain,  input-output  data  and  then  constructs  the  estimated  model  with 
the  Markov  parameters. 

Disturbance  rejection  algorithms.  As  mentioned  earlier,  one  of  the  possible  control 
schemes  for  closed-loop  separation  control  is  to  reduce  velocity  and  pressure  fluctuations  in  the 
separated  region.  This  control  scheme  is  generally  called  disturbance  rejection. 

Disturbance  rejection  controllers  have  been  widely  used  in  active  noise  control 
applications  (Kuo  and  Morgan  1996).  Recently,  researchers  have  started  to  apply  adaptive 
controllers  to  flow  control  problems.  For  example,  Cattafesta  et  al.  (1999)  used  an  adaptive 
system  to  suppress  the  disturbance  induced  by  the  flow  over  a  weapons-bay  cavity.  The 
advantages  of  using  adaptive  controllers  are  that  they  can  adapt  themselves  according  to  different 
flow  conditions  and  that  they  can  potentially  reduce  the  energy  cost  associated  with  the  flow 
control  problems.  Cattafesta  et  al.  (1997)  showed  that  the  control  of  cavity  flow  with  closed- 
loop  control  requires  one  order-of-magnitude  less  power  than  that  with  open  loop  control. 

Commonly  used  disturbance  rejection  algorithms  include  Filtered-X  LMS  (FXLMS), 
Filtered-U  LMS  (FULMS),  Filtered-X  RLS  (FXRLS)  and  Filtered-U  RLS  (FXRLS)  algorithms 
(Kuo  and  Morgan  1996).  Besides  these,  the  ARMARKOV  adaptive  control  algorithm  was  first 
introduced  by  Venugopal  and  Berstein  (1997)  and  further  developed  by  Sane  et  al.  (2001).  The 


17 


underling  model  structure  of  the  ARMARKOV  adaptive  control  algorithm  is  the  ARMARKOV 
representation,  which  is  an  extension  of  the  ARMA  representation  with  explicit  impulse  response 
(Markov)  parameters.  The  ARMARKOV  adaptive  control  algorithm  doesn’t  require  a  model  of 
the  control-to-reference  transfer  function  nor  does  it  require  a  model  of  the  transfer  function  from 
plant  disturbances  to  sensors  (Sane  et  al.  2001).  The  only  transfer  function  needed  is  the  control 
to  performance  transfer  function,  which  can  be  identified  simultaneously  using  the  recursive 
ARMARKOV/Toeplitz  system  identification  algorithm  described  in  the  previous  section. 

1.4  Objectives 

•  To  explore  suitable  linear  and  nonlinear  control  objectives  and  strategies  for 
closed-loop  control  of  separated  flows. 

•  To  implement  optimization  algorithms  and  system  identification/disturbance 
rejection  algorithms  for  closed-loop  control  of  separated  flow  on  a  wind  tunnel 
airfoil  model  (NACA  0025). 

•  To  analyze  performance,  adaptability,  costs,  and  limitations  of  closed-loop 
separation  control  algorithms. 

•  To  investigate  the  relevant  flow  physics  of  successful  feedback  control  strategies. 

1.5  Approach 

The  proposed  closed-loop  separation  control  includes  two  key  parts:  modeling  and 
control  strategies.  As  far  as  modeling  is  concerned,  two  types  of  approaches  can  be  implemented 
to  model  the  flow  characteristics:  1)  a  reduced-order  flow  model  based  on  the  Navier-Stokes 
equations,  2)  system  identification  techniques.  The  first  approach  is  widely  used  in 
computational  flow  control  simulations.  This  research  will  concentrate  on  experimental  studies 
by  using  system  identification  techniques  that  have  not  yet  been  applied  to  the  separation  control 
problem.  The  dynamical  systems  model  will  include  the  dynamics  of  actuators,  sensors,  and  the 
flow  system.  Then  the  disturbance  rejection  algorithm  is  used  to  suppress  flow  fluctuations  (e.g., 
measured  by  unsteady  pressure  transducers). 

One  the  other  hand,  for  the  non-model  based  optimization  algorithms,  no  system 
identification  is  needed.  The  possible  cost  functions  for  the  algorithm  are  summarized  as 
follows.  Since  the  suction  pressure  region  of  the  upper  surface  of  the  airfoil  is  primarily 
responsible  for  lift  generation  and  drag  reduction,  the  static  pressure  recovery  coefficient 
dCp/d(x/c)  over  the  upper  surface  of  the  airfoil  is  a  reasonable  candidate  as  a  cost  function  to 

maximize  for  feedback  separation  control.  Other  candidates  for  cost  functions  are  lift  and  drag 
or  combinations  of  these  (e.g.,  lift/drag  ratio).  The  benefit  of  using  lift/drag  is  that  L/D  is  a 
global  or  integrated  quantity  and  is  less  sensitive  to  sensor  location.  The  objectives  for  the 
controller  are  clear,  i.e.  to  minimize  drag  and  to  maximize  lift  or  the  ratio  of  lift/drag.  The 
experimental  setup  uses  a  lift/drag  balance  for  this  purpose. 

1.6  Outline  of  This  Dissertation 

A  theoretical  background  on  system  identification,  control,  and  optimization  algorithms 
will  be  discussed  in  Chapter  2.  Simulation  results  and  validation  experiments  of  the  algorithms 
will  then  be  presented  in  Chapter  3.  Chapter  4  describes  the  experimental  setup  and  techniques 
for  this  research.  Chapter  5  presents  experimental  results  and  discussion.  Summary  and  future 
work  will  be  presented  in  the  last  chapter. 
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Figure  1-1.  Separation  of  flow  over  an  airfoil. 
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Figure  1-2.  Types  of  velocity  profiles  as  a  function  of  pressure  gradient  (White  1991). 

Error!  Objects  cannot  be  created  from  editing  field  codes. 

Figure  1-3.  Lift  and  drag  coefficients  of  NACA  0025  airfoil  at  Re  =  100,000 . 
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Figure  1-5.  Characterization  of  possible  frequency  scales  in  separated  flow  (Mittal  et  al.  2005). 
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Figure  1-6.  Comparison  between  the  lumped  element  model  (-)  and  experimental  frequency 
response  (□)  measured  using  phase-locked  LDV  for  a  prototypical  synthetic  jet  (Gallas  et 
al.  2003). 
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2  Theoretical  Background 


This  chapter  presents  detailed  descriptions  and  derivations  of  the  algorithms  that  are  used 
in  this  research.  The  algorithms  include  optimization,  system  identification,  and  disturbance 
rejection  algorithms. 

2.1  Optimization  Algorithms 

Some  established  minimization  and  maximization  algorithms  are  summarized  by  Press  et 
al.  (1992).  The  downhill  simplex  algorithm  and  the  Powell’s  algorithm  do  not  require  derivative 
calculations,  which  makes  them  good  candidates  for  this  research  since  derivative  calculations 
are  problematic  for  (usually  noisy)  experimental  data.  Between  these  two  algorithms,  the 
downhill  simplex  algorithm  is  more  concise  and  self-contained.  The  so-called  extremum¬ 
seeking  algorithm  has  been  applied  to  a  flow  control  problem  by  Banaszuk  et  al.  (2003).  Thus, 
this  algorithm  is  also  summarized  here. 

2.1.1  Downhill  Simplex  Algorithm 

The  downhill  simplex  algorithm  is  implemented  to  minimize  an  objective  function  (e.g., 
drag-to-lift  ratio).  The  benefits  of  the  algorithm  are  its  simplicity,  applicability  to 
multidimensional  optimization  and  robust  performance.  The  algorithm  searches  downhill  in  a 
straightforward  fashion  that  makes  no  prior  assumptions  about  the  function.  The  downhill 
simplex  algorithm  requires  only  function  evaluations,  not  derivatives.  Since  it  does  not  make 
any  assumptions  about  the  function,  it  can  be  very  slow  sometimes.  However,  it  can  be  very 
robust  in  the  sense  that  it  guarantees  to  find  a  minimum  (at  least  a  local  minimum)  (Press  et  al. 
1992). 

A  simplex  is  the  geometrical  object  consisting,  in  N  dimensions,  of  N+l  points  (or 
vertices)  whereas  the  N+l  points  span  a  N -dimensional  vector  space  (Press  et  al.  1992).  For 
example,  in  two  dimensions,  a  simplex  is  a  triangle.  In  three  dimensions,  it  is  a  tetrahedron, 
although  not  necessarily  a  regular  tetrahedron.  The  downhill  simplex  algorithm  makes  use  of  the 
geometrical  concept  of  a  simplex  and  works  its  way  in  the  local  downhill  direction  until  it 
encounters  a  (at  least,  local)  minimum. 

The  key  steps  of  the  downhill  simplex  algorithm  are  summarized  as  follows: 

•  Evaluate  the  cost  function  at  chosen  initial  points.  Note  that  there  should  be  N+l 
initial  points,  defining  an  initial  simplex.  For  two  or  higher  dimensions,  the  initial 
points  should  not  be  linearly  dependent. 

•  Take  a  series  of  steps  to  move  in  the  downhill  direction.  As  an  example,  the  steps 
for  three-dimensional  search  are  illustrated  in  Figure  2-1.  In  the  figure, 
“Reflection”  means  that  the  algorithm  reflects  the  highest  (i.e.  worst)  point  about 
the  center  of  the  three  lower  (i.e.  better)  points  with  some  coefficients  to  the  other 
side  of  the  plane  and  then  evaluates  the  cost  function  at  the  reflected  point. 
“Expansion”  means  to  expand  further  along  the  reflection  direction  when  the 
“Reflection”  point  does  improve  (i.e.,  lower  the  cost  function).  “Contraction” 
means  to  move  the  highest  (i.e.  worst)  point  towards  the  plane  formed  by  the  three 
lower  (i.e.  better)  points,  thus  contracting  the  original  simplex.  To  summarize,  all 
the  necessary  steps  taken  here  are  to  move  the  worst  point  reference  to  the  plane 
formed  by  the  other  better  points  to  search  for  a  better  point. 
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•  Stop  when  some  termination  criteria  are  met.  For  example,  the  moving  distance 
is  smaller  than  some  tolerance  value. 

2.1.2  Extremum  Seeking  Algorithm 

Artiyur  and  Krstic  (2003)  present  the  theoretic  details  and  some  applications  of  the 
extremum-seeking  algorithm.  Simulations  using  the  algorithm  can  be  done  following  the  block 
diagram  in  Figure  2-2.  The  simple  proof  that  this  algorithm  will  drive  f(9)  to  its  extremum  is 
summarized  below. 

First,  assume  that  f(0)  has  a  minimum  f  *  and  can  be  approximated  as  the  following  form: 


•  \2 


2! 

where  0*  is  the  optimal  input  and  f  is  the  local  curvature  of  the  cost  function  f(0)  near  0* . 
Since  it  is  assumed  that  f(0)  has  a  minimum,  f  should  be  larger  than  zero  for  this  case. 
Next,  define  the  estimated  error: 

0=0*  -0 


(4) 


where  9  is  the  estimated  optimal  input. 

From  Figure  2-2, 

0=0+asin(wt)=0*  -0+asin(wt) 
Substituting  equation  (6)  into  equation  (4)  results  in 

y=f(0)=f  *  +  [asin(wt)-0  J 

„  ,  v  l-cos(2wt) 

Expand  equation  (7)  and  apply  sin  (wt)= -  to  obtain 


.  fa  .  2 

y=f  + 


2  2  -  f  02 

sin  (wt)-f  a0sin(wt)+ - 


,  f  a2  f  a2  f  02 

=f  + - cos(2wt)-f  a0sin(wt)+ - 


(5) 

(6) 
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From  Figure  2-2,  this  signal  y  will  pass  through  a  high  pass  filter 


s+wt 


.  Let  0<wh<w , 


then  all  the  DC  components  in  equation  (8)  will  be  removed  while  the  oscillatory  terms  remain. 

(9) 


f'a2 

q» - cos(2wt)-f  a0sin(wt) 

4 


Next,  77  is  multiplied  by  asin(wt)  to  give 

f'a2 


e= - cos(2wt)sin(wt)-f  a0sin(wt) 

4 


(10) 


Using  the  trigonometric  identities  sin2(wt)= 


l-cos(2wt) 


and 


...  .  sin(3wt)-sin(wt) 

cos(2wt)sin(wt)= - results  in 
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(11) 


f  a  sin(3wt)-sin(wt)  ~l-cos(2wt) 
£= - 1  at) - 


4  2 

f  a§  f  a2 


f  a 


f  a0 
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sin(wt) - sin(3wt)+ - cos(2wt) 

8  2 


From  Figure  2-2,  this  signal  e  passes  through  a  low  pass  filter.  Let  0<w,<w  ,  then  all  the 
high  frequency  terms  will  be  removed  and  only  the  DC  term  remains. 
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This  signal  then  passes  through  an  integrator 


This  gives 


e=— i 
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From  equation  (5),  assuming  f  *  is  fixed,  then 

0=-0 

From  equations  (14)  and  (15),  one  can  obtain  the  first-order  differential  equation 


with  the  solution 
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(15) 
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(17) 


Since  f  is  assumed  to  be  positive  and  a  and  k  are  positive  constants,  the  estimated  error 
0  will  exponentially  decay  to  zero. 

2.2  System  Identification  Algorithms 

System  identification  (ID)  uses  measured  signals  (i.e.,  inputs  and  outputs  of  the  system) 
to  identify  (or  estimate)  the  unknown  parameters  of  an  assumed  dynamical  systems  model.  It 
thus  provides  the  necessary  system  information  for  control  algorithms.  System  identification 
algorithms  can  be  divided  into  two  categories:  offline  (or  batch)  and  online  (or  recursive). 
Offline  algorithms  first  digitize  a  data  record  and  then  try  to  estimate  the  system  using  these  data 
offline,  usually  via  a  least  squares  method.  Conversely,  online  algorithms  identify  systems 
recursively  while  acquiring  data  in  real-time.  Online  system  identification  is  also  known  as 
adaptive  filtering. 

In  this  research,  three  system  ID  algorithms  will  be  investigated:  ARMARKOV/LS, 
ARMARKOV/LS/ERA  and  recursive  ARMARKOV/Toeplitz  algorithms.  They  are  all  based  on 
the  ARMARKOV  representation,  which  explicitly  contains  Markov  parameters  (i.e.,  pulse 
response)  of  the  system.  The  well  known  ARMA  representation  contains  only  one  Markov 
parameter  and  is  a  special  case  of  the  ARMARKOV  representation.  The  main  advantage  of 
these  algorithms  is  their  robustness  with  respect  to  low  signal-to-noise  ratios  (Akers  and 
Bernstein  1997  A,  B).  The  ARMARKOV/LS  algorithm  is  an  offline  algorithm  and  implements 
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an  overparameterized  realization  of  the  system.  The  ARMARKOV/LS/ERA  algorithm  uses  the 
same  procedures  to  identify  the  system  parameters  as  the  ARMARKOV/LS  algorithm,  but 
implements  a  minimal  realization  of  the  system.  The  recursive  ARMARKOV/Toeplitz 
algorithm  is  an  online  algorithm.  The  advantage  of  using  an  online  algorithm  is  that  it  can  adapt 
to  the  changing  system. 

2.2.1  ARMARKOV/LS  Algorithm 

Consider  the  discrete-time  finite-dimensional  linear  time-invariant  system: 

x(k+ 1  )=Ax(k)+Bu(k) 
y(k)=Cx(k)+Du(k) 


where Ae  Rnxn,Be  Rnx',Ce  Rlxn,De  R|XI ,  and  i  and  1  are  the  number  of  inputs  and  outputs, 
respectively,  of  the  system.  For  a  single-input/single-output  (SISO)  system,  i=l=  1 .  The 
algorithm  is  derived  below  for  a  SISO  system. 

The  Markov  parameters  Hj  are  defined  by 


Hj  =  D  j=-l 

H.  =  CAJB  j>0 

J  J 


(19) 


Next,  define  the  ARMARKOV  regressor  vector  d>  (k)e  R 


2n+n 


oM(k)^ 


y(k-p) 

y(k-p-n+l) 

u(k) 


(20) 


u(k-p-n+l) 


where  n  is  the  order  of  the  system  and  p  is  the  number  of  Markov  parameters.  Here,  y  and  u 
denote  measured  input  and  output  of  the  system  described  in  equation  (18),  respectively. 

Next  define  the  estimated  output  of  the  system 

y(k)=WOM  (21) 


where  the  ARMARKOV  weight  matrix  W  is  defined  by 


W=[-A„  H,LHm.2  BJ 


(22) 


(23) 


The  expression  of  the  weight  matrix  Wp  is  then  determined  to  minimize  the  output  error 

cost  function  defined  below. 

First,  define  the  output  error 

£(k) =  y(k)-y(k)  (24) 

and  the  output  mean  squared  error  cost  function 
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(25) 


1  N  1  1  N  1 

J  =  -  Y -  e2  (k)= — Y -  eT  (k)e(k) 

NtT2  N£f2 

where  N  is  the  number  of  measurements. 

Substituting  equations  (20),  (21),  (22)  and  (24)  into  equation  (25)  results  in 

J=-^X^(y(k)-W,»(k))T(y(k)-WO(k))  (26) 

IN  k=|  Z 

J  =  -^rS^(yT  (k)-<I>T  tk)WT)(y(k)-Wd>(k))  (27) 

IN  k=,  Z 

»,=ij£|<yT(|Oy(k)-®T(k)wTy(io 

^  k=i  ^  (28) 

-yT(k)W<t(k)+<DT(k)WTW(k)<D(k)) 

Because  <I)T  (k)  WTy(k)  and  yT(k)W<I>(k)  are  transposes  of  each  other  and  are  also 
scalars,  they  are  equal  to  each  other.  So, 


j=~X^(yT(k)y(k)-2<I>T(k)*Ty(k)+<,>T(k)wTw{k)a>) 

IN  k=,  Z 


(29) 


To  find  the  W  to  minimize  the  output  error  cost  function  defined  in  equation  (25),  we  set 
the  partial  derivative  of  J  with  respect  to  W  equal  to  zero.  So, 

9J  3rJ 


aw  awT  0 

From  matrix  calculus,  we  derive  each  term  of  J  in  equation  (29)  first, 


^(yT<k)y(k)M 


^f(2<I)T(k)  WTy(k))=2((0T  (k)y(k))T  )T  =2<kT  (k)y(k) 

3(0T(k)WTWD(k))d(wTW> 

3(WTW)  dWT 

=  [<DT(k)<D(k)(2WT)J 
=  2WOT(k)<F(k) 


aw 

(k)  WTWO(k))  = 


(30) 

(31) 

(32) 


(33) 


Thus, 


aTj 


j^  =  0  =  ^|;i[-2<I)T(k)y(k)+2W®T(k)<l>(k)] 


aw 


i-|;w4>T(k)<i)(k)=i|;®T(k)y(k) 

IN  k=i  IN  w- 1 


(34) 

(35) 
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Finally,  the  expression  of  the  weight  matrix  W  to  minimize  the  output  error  cost  function 
is  given  by 

-i-i 


W= 


-J-J>T(k)y(k) 
N  k=l 


-ir^0T(k)0(k) 

.  ^  k=l 


(36) 


After  extracting  the  coefficients  Am,B^  and  Hj  from  W  via  equation  (22),  we 

can  obtain  the  system  transfer  function  of  the  ARMARKOV  representation,  which  is  defined  as 
follows 


G,(z)= 


H,lz^-|+L+H(,,2zn+P>uzn-|+L+P^n 
zM+n‘‘  +a^,zn'‘  +L+aM 


(37) 


This  is  called  the  ARMARKOV/LS  identification  algorithm  and  this  algorithm  assumes 
the  numerator  has  the  same  order  of  the  denominator  for  simplicity.  For  the  systems  whose 
numerators  and  denominators  do  not  have  the  same  order,  some  parameters  described  in  equation 
(37)  will  be  identified  to  be  approximately  zero. 

The  well-known  ARMA  representation  only  has  one  explicit  Markov  parameter  and  it  is  a 
special  form  of  the  ARMARKOV  representation  with  p=l  in  (37) 


G2(z)= 


H-|Zn+P|  1zn~l+L+P,n 
zn+a,  ,zn  l+L+al  n 


(38) 


For  system  identification  problems,  the  order  of  the  system  is  usually  not  known  in 
advance,  so  we  adjust  n  and  //  to  improve  the  performance  of  the  system  identification 
algorithm. 

2.2.2  ARMARKOV/LS/ERA  Algorithm 

The  ARMARKOV/LS/ERA  algorithm  obtains  a  minimal  realization  of  the  transfer 
function  of  the  system  from  the  Markov  parameters.  It  uses  the  same  algorithm  as  the 

ARMARKOV/LS  algorithm  to  obtain  the  weight  matrix  W  by  equation  (36).  Then  the  Markov 
parameters  H(  can  be  extracted  from  equation  (36)  by  using  equation  (22).  Next,  define  the 


Markov  block  Hankel  matrix  for  a  SISO  system: 


H,sj  = 


'Hi  -  HJ+s" 


H  . . .  H 
V  j+r  j+r+s/ 


(39) 


where  r,s  are  any  positive  integers.  In  this  research,  r  is  set  to  be  equal  to  s  for  convenience. 

Then,  we  apply  the  singular  value  decomposition  as  describe  in  Akers  and  Bernstein 
(1997 A) 


H„0=PSr,QT 


(40) 


where  P1P=QIQ=I  and  Srs  =  diagonal  matrix  of  singular  values. 

From  the  Eigenvalue  Realization  Algorithm  (ERA),  (Juang  1994,  pp.  133-137) 
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-1/2 


where  E.  = 


IX I 


A=Sr;,/2PTHrAlQSr, 
B=Srs,/2QTEs 
C=EtPS  1/2 

r  r,s 

D=H  , 

.  The  ERA  also  requires  r,s  >  n-1 . 


This  arrives  at  the  minimal  state  space  realization  of  the  system 

G3(z)=C(zA-I)  '  B+D 


(41) 


(42) 


This  is  called  the  ARMARKOV/LS/ERA  algorithm.  It  is  a  minimal  realization  because 
the  system  order  can  be  chosen  as  a  minimal  value  when  using  the  singular  value  decomposition 
in  equation  (40).  However,  an  important  drawback  of  this  algorithm  is  that  the  singular  value 
decomposition  in  equation  (2.23)  is  very  computational  intensive. 

Theoretically,  the  rank  of  the  Srs  matrix  should  be  the  rank  of  the  system.  However,  in 

practical  applications,  the  singular  value  decomposition  will  return  more  singular  values  than  the 
system  order  due  to  measurement  noise,  and  so  the  extra  singular  values  should  be  small.  So, 
only  the  largest  n  singular  values  obtained  by  the  singular  value  decomposition  will  be  used. 


2.2.3  Recursive  ARMARKOV/Toeplitz  Algorithm 

First,  define  the  ARMARKOV  regressor  vector  d>M  (k)  e  r2,1+2p-2+m  ; 

y(k-p) 


<I>M(k)  = 


y(k-p-p-n+2) 

u(k) 


u(k-p-p-n+2) 


(43) 


where  n  is  the  order  of  a  system,  p  is  the  number  of  Markov  parameters,  and  a  new  parameter 
p  determines  the  averaging  window  of  input-output  data  that  appears  in  the  above  regressor 
vector. 

It  follows  that 


y(k)=W<DM  (44) 

where  the  ARMARKOV/Toeplitz  weight  matrix  We  rpx<2"+2p-2+p>  js  the  block-Toeplitz  matrix 
defined  by 

'A  0  -  0  H,  -  h,.2  bm  0-0 

o  :  0  ■•.  : 

:  o  ;  o 

•••  o  H, 


WM  — 


o  •••  o  a  0 


HM-2  B, 


(45) 
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and  Am  =  [apl  •••a(1>n]I  R  *" ,  B(1  =  (Ppl  •••?„,„]€  R  ,  and  Hj  are  the  Markov  parameters. 
As  before,  define  the  output  error  e(k)  and  the  output  error  cost  function  J  (k) 

e(k)  =  Y(k)-Y(k) 


J(k)  =  ^eT(k)e(k) 

Next,  the  gradient  of  J  (k)  with  respect  to  W(k)  can  be  calculated  by 

^P-=-uore(k)d)T(k)i 
aw(k)  L  J 

where  °  denotes  the  Hadamard  product  (i.e.  element-wise  matrix  product)  and  Ue  R 
is  defined  by 

0 


(46) 

(47) 

(48) 

px(2n+2p-2+n) 


U  = 


I,-  0 

0 


0 


®  ^lx(ji+n)  ® 


0 


0 


0 


0  I,»0 


0  I 


lx(p+n) 


(49) 


Finally,  the  recursive  update  law  for  the  weight  matrix  W  is  given  by 

W(k+l)  =  W(k)-ri(k)-aj(k) 


(50) 


aw(k) 

In  equation  (50),  q(k)  is  the  adaptive  step  size.  The  optimal  adaptive  step  size  nopt(k)  is 


defined  as 


-ooE 


aj(k) 


aw(k) 


(51) 


where  denotes  the  spectral  norm. 

The  computationally  efficient  step  size  qeff  (k)  (namely,  it  is  more  computational  efficient 
since  it  only  needs  to  calculate  the  normal  ARMARKOV  regressor  vector  (k) )  is  defined  as 

1 


n,„00  = 


(52) 


In  order  to  assure  convergence,  q(k)  should  satisfy  r|(k)=ar|opt  (k)  or  q(k)=aqcff  (k) , 
where  a  e  (0,2) . 

After  W  matrix  is  obtained  by  (50),  we  can  extract  the  coefficients  Ap,Bp  and  Hj  from 

(22).  Then,  we  can  obtain  the  system  transfer  function  of  the  ARMARKOV  representation  form, 
which  is  defined  in  equation  (37).  Since  the  A(i,Bp  and  Hj  coefficients  are  updated  every 

iteration,  this  algorithm  is  called  as  the  recursive  ARMARKOV/Toeplitz  algorithm. 
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2.3  Adaptive  Disturbance  Rejection  Algorithms 

Disturbance  rejection  controllers  have  been  widely  used  in  active  noise  control 
applications  (Kuo  and  Morgan  1996).  A  block  diagram  of  a  standard  disturbance  rejection 
problems  is  shown  in  Figure  2-3,  where  w  is  the  disturbance,  u  is  the  control  signal,  y  is  the 

reference  signal,  z  is  the  performance  signal  and  Gc  is  the  disturbance  rejection  controller.  The 
goal  for  the  controller  is  to  generate  a  control  signal  u  to  minimize  some  cost  function  of  the 
performance  signal.  The  four  transfer  matrices,  namely,  the  primary  path  Gzw ,  the  secondary 

path  Gzu,  the  reference  path  G^  and  the  feedback  path  G^,  are  standard  terminology  in  the 

noise  control  literature  (Kuo  and  Morgan  1996).  The  feedforward-type  disturbance  rejection 
algorithms,  such  as  FXLMS  and  FXRLS,  assume  that  G^O  (no  feedback  path)  and  G^I. 

On  the  other  hand,  the  ARMARKOV  disturbance  rejection  algorithm  does  not  make  these 
assumptions.  All  the  disturbance  rejection  algorithms  require  identifying  the  secondary  path  Gzu 

by  online  or  offline  system  identification  methods. 

2.3.1  ARMARKOV  Disturbance  Rejection  Algorithm 
Consider  the  linear  discrete  time  two-input/two-output  system  (shown  in  Figure  2-3)  given 
by 

z(k)=Gzww(k)+Gzuu(k)  (53) 

y(k)=Gyww(k)+Gyuu(k)  (54) 

where  the  disturbance  w(k) ,  the  control  u(k) ,  the  reference  y(k )  and  the  performance  z(k)  are 

in  Rmw  ,  Rm”,  R ,y  and  R1,  respectively,  and  m  and  1  denote  the  number  of  inputs  and  outputs, 
respectively.  The  system  transfer  matrices  Gzw  (primary  path),  Gzu  (secondary  path),  Gyw 

(reference  path)  and  G^  (control  path)  are  in  Rl,xmw ,  Rl,xn\  R|yxn’"  an(j  R'>xm“  ^  respectively. 
The  objective  of  the  active  noise  or  vibration  control  problems  is  to  determine  a  controller 
Gc  e  Rm'xl’  that  produces  a  control  signal  u(k)=Gcy(k)  such  that  the  performance  measure  z(k) 
is  minimized  (Sane  et  al.  2001).  A  measurement  of  z(k)  is  used  to  adapt  Gc . 

The  ARMARKOV  form  of  (53)  -  (54)  is 

z(k)  =  -ajZ(k-p-j-l  )+£  HZWJ.2w(k-j+ 1  )+£  BZWJw(k-p-j- 1 ) 

j;'  j='n  H  (55) 

+X  H*.j-2U(k-j+ 1  )+X  Bz«jU(k-F-j+ 1 ) 

j=l  j=! 


y(k)=]T -ajy(k-p-j-l  )+£  Hywj.2w(k-j+l  )+£  Bywjw(k-p-j-l ) 

j=l  j=l  j=l 

+S  H^u(k-j+l)+ £  B„jU(k-H+ 1 ) 

j=l  j=l 


(56) 


where  R,  B 


zwj 


R 


Lxm„, 


BZu,j,HZUJ 


G  R 


Lxm., 


r>  o  r- 

Dywj’°ywJ 


R 


lvxmu. 


^yuj’^yuj  ' 


p|)xni. 


is  the  order  of  the  system,  and  is  the  number  of  the  Markov  parameters. 


n 
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Then,  we  define  the  extended  performance  vector  Z(k),  the  extended  measurement 
vector  Y(k),  and  the  extended  control  vector  U(k)  as 


Z(k)  =  [z(k)  L  z(k-p+l)]T  (57) 

Y(k)  =  [  y(k)  L  y(k-p+ 1  )]T  (58) 

U(k)  =  [u(k)  L  u(k-pc+l)]T  (59) 

where  p  is  an  averaging  or  windowing  parameter  and  pc=(p+n+p-l) . 

The  ARMARKOV  regressor  vectors  <DZW  and  are  defined  by 

<DZW  =  [  z(k-p)  •  •  •  z(k-p-p-n+2)  w(k)  •  •  •  w(k-p-p-n+2)]T  (60) 

and 

<Dyw  -  [  y(k-p)  ■  •  •  y(k-p-p-n+2)  w(k)  •  •  •  w(k-p-p-n+2)] r  (61) 

Then  (55)  and  (56)  can  be  written  as 

Z(k)=WzwOzw  +BzuU(k)  (62) 

Y(k)=WywOyw+ByuU(k)  (63) 


where  Wzw ,  Bzu ,  Wzw  and  Byu  are  the  ARMARKOV  weight  matrices.  Only  Bzu  will  be  used  in 
the  control  algorithm  (shown  later),  and  it  will  be  obtained  using  the  ARMARKOV/Toeplitz 
system  identification  algorithm.  The  ARMARKOV  control  matrix  Bzu  is  given  by 


B  = 

zu 


Hz,,  ^zu,n-2  BzuJ  Bzun  0[  xrnlJ 

0, 


l,xm, 


0 


l.xrn.. 


Hzu,.i  HZU|1.2  Bzu ,  •••  Bzun 


(64) 


where  0,  xm  is  the  zero  matrix. 

Next,  the  ARMARKOV  adaptive  disturbance  rejection  algorithm  is  derived.  The  control 
signal  u(k)  is  given  by 


u(k)=£  -ac  jU(k-pc  -j+ 1 )+ £  HCj2y(k-j+ 1  )+£  Bcjy(k-pc  -j+ 1 ) 
j=i  j=i  j=i 

where  a  . e  R  and  H  ,,B cie  Rm"x1’ . 

Similarly,  the  delayed  versions  are 

u(k- 1  )=£  -acJu(k-pc  -j)+£  HcJ.,  y(k-j)+£  Bcjy(k-pc  -j) 
i=i  j=i  j=i 


(65) 


(66) 


31 


(67) 


J'c 

u(k-pc  + 1  )=X  ■acju(k’^c  -J-Pc  +2) 


J=» 

He 


+ZHcj-iy(k-j-Pc+2)+ZBcjy(k-Pc-j-pc+2) 

j=i  j=i 

Substituting  all  these  equations  in  (59)  and  reordering  gives 

U(k)=|]Li0(k-i+l)Ri<I>uy(k) 


where 


9(k)i[-a,,(k)I„_  aCB<  (k)Im>  Hc„(k)  ■Ht,_.i(k)  Bc,,(k)- B,„(k)] 

°uy  (k)  =  [ u(k-pc )  •  •  •  u(k-pc  -nc  -pc  +2)  y(k)  •  •  •  y(k-pc  -nc  -pc  + 1  )]T 


L;  = 


B(i-I)m11xmn 


0 


and 


R,  = 


(pc  •ijm.xm, 


Oq,x(i-l)m„  ^q,xq,  ^q,x(pc-i)mu  ^q,x(i-l)ly  ^q,xq2  ^q,  x(pc -i)ly 
Oq2x(j-l)m„  ^q2x4.  ®q2x(Pc>>m0  ®q2x(i-l)l  ^q2xq2  ^q2x(pc-i)l 


(68) 

(69) 

(70) 

(71) 

(72) 


with  q,  =  ncmu  and  q2  =  (nc+pc-l)ly . 

Thus  from  (62)  and  (68),  we  obtain 

Z(k)  =  WlwOtw(k)+BIu^L,9(k-i+l)R1<J>uy(k)=Wrw(I>rw(k)+B!uU(k) 


(73) 


i=l 


Next,  evaluate  the  performance  of  the  current  value  of  0(k)  based  upon  the  behavior  of 
the  system  during  the  previous  p  steps  to  result  in  the  definition  of  the  estimated  performance 

Z(k)  by 

Z(k)  ^  Wzwd)zw(k)+BzuXL10(k)R,(I)uy(k)  (74) 

i=l 

Substituting  (73)  into  (74),  we  obtain  the  estimated  performance  in  terms  of  known  and 
measured  variables 


Z(k)  =  Z(k)-Bz 


U(k)-XL,0(k)R,d)uy(k) 


i=l 


Using  (75),  we  define  the  estimated  performance  cost  function 

J(k)=^ZT(k)Z(k) 


(75) 


(76) 
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The  purpose  for  the  ARMARKOV  adaptive  controller  is  to  obtain  the  controller 
parameters  0(k)  such  that  the  performance  cost  function  J(k)  is  minimized.  Using  matrix 
derivative  formulae,  the  gradient  of  J(k)  with  respect  to  0(k)  is  given  by 


9J(k) 

30(k) 


=  ZLlB;,Z(k)4.;<k)R^ 


The  gradient  is  used  in  the  update  law 


0(k+l)=0(k)-ti(k) 


3J(k) 

30(k) 


(77) 


(78) 


where  r|(k)  is  the  adaptive  step  size.  An  implementable  adaptive  step  size  rjmp{k)  is  used 


(79) 


where  and  denote  the  Frobenius  norm  and  the  spectral  norm  (Golub  and  Van  Loan 
1996),  respectively. 

The  steps  involved  in  implementing  the  ARMARKOV  adaptive  disturbance  rejection 
algorithm  are  summarized  as  follows: 

•  Obtain  the  matrix  Bzu  (eqn.  (64))  by  using  the  recursive  ARMARKOV/Toeplitz  system 
identification  algorithm  (eqn.  (45))  or  the  offline  ARMARKOV/LS  (eqn.  (36)). 

•  Calculate  the  control  signal  u(k)  from  the  controller  parameter  matrix  0(k)  and  the  vector 
4>uy(k)  (eqn.  (68)). 

•  Use  the  signals  u(k),  z(k)  and  y(k)  to  update  the  vectors  Z(k)  (eqn. (75))  and  Ouy(k) 

(eqn.  (70)). 

•  Calculate  the  gradient  —  (eqn.  (77)). 

O0(k) 

•  Calculate  the  implementable  adaptive  step  size  r|imp(k)  (eqn.  (79)). 

•  Update  the  controller  parameter  matrix  0(k)  ->  0(k+l)  (eqn  (78)). 
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Figure  2-1 .  Flow  chart  of  downhill  simplex  algorithm. 
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Figure  2-2.  Block  diagram  for  the  extremum  seeking  control. 


Figure  2-3.  Block  diagram  of  disturbance  rejection  control. 
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3  Simulation  and  Validation  experiments 


Before  the  algorithms  are  used  for  closed-loop  separation  control  in  the  wind  tunnel 
experiments,  they  are  tested  by  using  Matlab/Simulink  simulations  or  validation  experiments. 
The  purpose  of  this  chapter  is  to  ensure  that  the  algorithms  work  as  desired. 

3.1  Optimization  Simulations 

3.1.1  Downhill  Simplex  Simulation  Results 

The  downhill  simplex  algorithm  is  programmed  in  Matlab.  The  performance  of  the 
algorithm  is  illustrated  by  a  1 -dimensional  and  a  2-dimensional  simulation  cases.  This  algorithm 
can  be  easily  extended  to  higher  dimensions. 

In  the  1 -dimensional  case,  the  cost  function  f(x)  is  chosen  to  be  an  8th  order  polynomial 

function  of  x ,  which  has  a  local  minimum  at  x=14.2  and  a  global  minimum  at  x=67.3  as  shown 
in  Figure  3-1.  Two  initial  conditions  are  selected  to  demonstrate  that  this  algorithm  can  be 
“trapped”  by  a  local  minimum.  The  first  initial  condition  is  at  about  x=45.  One  should  notice 
that  for  this  1 -dimensional  problem,  there  should  be  two  independent  points  (a  simplex)  as  the 
initial  condition.  As  shown  in  Figure  3-1,  the  downhill  simplex  algorithm  crawls  down  to  the 
global  minimum  (red  trace).  On  the  other  hand,  the  second  initial  condition  is  at  about  x=30 , 
which  leads  the  algorithm  to  the  local  minimum  (blue  trace).  This  is  dictated  by  the  inherent 
“downhill”  nature  of  the  algorithm. 

Another  example  is  to  demonstrate  how  the  downhill  simplex  algorithm  works  in  two- 
dimensional  space.  The  cost  function  is  obtained  in  Matlab  by  the  “peaks”  command.  The 
formula  for  the  cost  function  is  as  follows: 

Z  =  3  [( 1  -x)2e',!  -(y+ 1  )2  ]  - 1 0(  -x3  -y5  )e*2  'yJ  -  ^  e^x+,)V  (80) 

This  function  has  two  local  minima  and  a  global  minimum  as  shown  in  Figure  3-2. 
Similar  as  the  1-dimensoinal  case,  the  optimization  algorithm  converges  to  either  a  local 
minimum  (blue  trace)  or  the  global  minimum  (red  trace)  depending  on  the  initial  condition. 
Although  each  iteration  of  the  algorithm  requires  several  steps  (Chapter  2),  it  only  takes  9 
iterations  to  find  the  global  minimum.  This  result  is  encouraging  and  suggests  that  it  can  be  fast 
for  some  cases.  One  can  also  adjust  the  termination  tolerance  to  control  the  time  consumption. 
On  the  other  hand,  the  time  consumption  of  the  separation  control  experiments  is  also  dependent 
on  other  factors,  such  as  data  acquisition.  This  will  be  discussed  further  in  Chapter  5. 

3.1.2  Extremum  Seeking  Simulation  Results 

The  extremum  seeking  algorithm  is  implemented  in  Simulink.  Figure  3-3  shows  the 
simulation  block  diagram  for  the  extremum  seeking  control.  In  the  simulation,  the  algorithm 
seeks  a  maximum  instead  of  a  minimum.  One  can  easily  modify  the  program  to  search  for  a 
minimum  by  adding  a  negative  sign  to  the  cost  function.  Two  numerical  models  are  tested.  The 

first  model  is  a  quadratic  function  f=f *-(0-0* )  ,  which  has  a  single  maximum  f*  at  0=0’ as 

shown  in  Figure  3-4.  In  this  case,  f*  is  set  to  10  and  0*  is  set  to  5.  The  second  model  is  a 
double  hump  model,  which  is  fitted  by  a  8th  order  polynomial  function,  which  is  the  same  as  the 
model  shown  in  Figure  3-1  with  a  opposite  sign.  It  has  a  local  maximum  and  a  global  maximum 
as  shown  in  Figure  3-5. 
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The  ARMARKOV  disturbance  rejection  algorithm  requires  the  following  three 
parameters:  the  order  of  the  controller  nc ,  the  number  of  Markov  parameters  of  the  controller  pc 
and  the  adaptive  step  size  constant  y  that  controls  the  convergence  rate  of  the  controller. 

The  controller  uses  the  system  information  identified  by  the  ARMARKOV/Toeplitz 
system  ID  with  the  parameters  of  p=n=2,  p=10  and  SNR=20  dB  (shown  in  Figure  3-21).  All 
cases  of  the  controller  design  use  the  same  identified  system. 

Band-limited  white  noise  with  frequency  of  0  -  150  Hz  and  variance  of  0.09  is  used  to 
excite  the  disturbance  piezoceramics.  Figure  3-26  and  Figure  3-27  show  the  time  data  of  the 
performance  and  control  signals  with  the  parameters  of  nc=2 ,  pc=20  and  T  =  1 .  The  system  ID 

is  off  for  the  whole  period  and  the  controller  is  off  initially  and  turned  on  at  t  =  20  sec.  In  Figure 
3-28,  the  power  spectra  of  the  performance  signal  with  control  off  and  on  are  compared.  The 
power  spectra  were  calculated  by  using  the  time  data  of  20-second  duration  with  NFFT=1024, 
50%  overlap  and  a  hanning  window.  The  “control  on”  case  is  taken  after  the  controller  is  turned 

on  for  30  seconds.  The  performance  of  suppression  is  calculated  by  10glog10  (o20ff/c2on )  and  this 

case  gives  1 1.7  dB  suppression.  Interestingly,  in  Figure  3-28,  the  power  around  70  Hz  and  120 
Hz  of  the  “control  on”  case  is  higher  than  that  of  the  “control  off’  case.  This  is  generally  defined 
as  spillover  (Hong  and  Bernstein  1998).  Hong  and  Bernstein  (1998)  used  the  Bode  integral 
constraint  to  analyze  the  spillover  problem  and  concluded  that  the  spillover  is  inevitable  if  the 
reference  and  performance  signals  are  collocated  or  the  disturbance  and  control  actuators  are 
collocated.  For  this  vibration  control  test,  the  reference  and  performance  signals  are  collocated, 
thus  the  spillover  is  unavoidable. 

Figure  3-29,  Figure  3-30  and  Figure  3-32  show  the  performance  of  the  disturbance 
rejection  algorithm  with  varying  nc ,  pc  and  T ,  respectively,  while  other  parameters  are  held 
constant.  From 

Table  3-2  and  Table  3-3,  it  is  interesting  to  find  that  there  is  not  much  difference  of  the 
suppression  performance  for  varying  nc  and  pc .  However,  the  step  size  parameter  T  does  play 
a  significant  role.  Larger  T  gives  much  faster  convergence  and  better  performance.  However, 
if  T  is  too  large,  it  is  possible  for  the  controller  to  become  unstable.  This  tradeoff  should  be 
kept  in  mind  when  choosing  T  . 

The  disturbance  rejection  controller  can  also  be  run  at  the  “ED  and  control”  mode.  In  this 
mode,  the  band-limited  white  noise  with  frequency  of  0  -  150  Hz  and  variance  of  0.09  is  used  to 
excite  the  disturbance  piezoceramics.  Meanwhile,  the  band-limited  white  noise  with  frequency 
of  0  -  150  Hz  and  variance  of  0.01  is  added  to  the  controller  output.  The  control  signal  is  shown 
in  Figure  3-33  and  the  performance  signal  is  shown  in  Figure  3-34.  Comparing  with  Figure 
3-27,  the  control  signal  of  the  “ID  and  control”  case  is  significantly  larger  at  the  beginning 
because  of  the  additive  signal  uID  and  the  evolution  of  the  controller  output  is  “buried”  under  it. 

Figure  3-35  compares  the  power  spectra  of  the  performance  signal  of  the  two  different 
modes.  It  is  surprising  that  the  “ID  and  control”  mode  results  in  lower  power  around  the  natural 
frequency  of  the  beam.  This  is  hard  to  see  in  Figure  3-34  because  the  “ID,  then  control”  mode 
seems  much  better.  However,  it  is  not  surprising  that  the  “ID  and  control”  mode  results  in  higher 
power  at  other  frequencies  than  the  natural  frequency  because  of  the  additive  signal  uID.  In 
addition, 

Table  3-5  shows  that  the  “ID,  then  control”  mode  gives  better  suppression  performance 
of  the  overall  power. 
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Unfortunately,  for  this  setup,  it  is  not  feasible  to  test  the  adaptability  of  the  two  modes. 
However,  this  will  be  done  in  the  wind  tunnel  experiments. 

As  a  summary,  the  computational  tests  are  conducted  first  to  determine  how  the 
parameters  affect  the  computational  complexity  of  the  system  ID  and  control  algorithms.  It  is 
shown  that  the  averaging  window  number  p  has  much  more  significant  impact  on  the 
computational  intensity  than  the  other  two  parameters  for  the  system  ID  algorithm.  The 
dependence  of  the  computational  complexity  vs.  n  is  approximately  twice  of  that  for  p . 
Similarly  it  is  found  that  the  turnaround  time  of  the  control  algorithm  is  approximately  linearly 
proportional  to  nc  and  pc ,  while  the  slope  for  nc  is  approximately  twice  of  that  for  pc . 

The  ARMARKOV/Toeplitz  system  ID  algorithm  successfully  identifies  the  system 
(control  model)  and  results  in  very  good  frequency  response  approximations.  A  significant 
improvement  of  the  performance  of  the  ARMARKOV  system  ID  over  the  ARMA  (when  p=l ) 
system  ID  is  found.  However,  too  many  Markov  parameters  of  the  ARMARKOV  system  ID 
may  be  detrimental  to  the  performance.  Higher  SNR  improves  the  performance,  thus  when  the 
system  ID  is  conducted  with  unknown  noise,  the  input  signal  should  be  chosen  as  large  as 
possible  within  the  maximum  allowable  level. 

The  order  of  the  controller  nc  and  the  number  of  Markov  parameters  pc  do  not  play 

significant  roles  on  the  performance  of  the  ARMARKOV  controller  for  this  vibration  control 
test.  However,  this  conclusion  may  vary  with  different  systems  and  remains  to  be  investigated. 
The  step  size  constant  T  significantly  affects  the  convergence  rate  of  the  controller.  T  should 
be  chosen  as  large  as  possible  before  it  makes  the  controller  unstable.  The  spillover  effect  is 
identified  in  this  vibration  control  test.  This  effect  is  unavoidable  because  the  “reference  signal” 
and  “performance  signal”  are  collocated  (Hong  and  Bernstein  1998). 

summarizes  all  the  parameters  that  are  used  in  the  simulations.  Recall  that  the  detailed 
derivation  of  algorithm  is  given  in  Chapter  2.  Note  that  the  parameters  a  and  w  are  the  main 
factors  that  affect  the  convergence  rate  and  stability.  Thus,  they  are  varied  in  the  simulations  to 
understand  how  they  affect  the  performance  of  the  algorithm. 

Figure  3-6  demonstrates  how  a  affects  the  convergence  rate  while  w  is  fixed  to  be  50 
Hz.  Clearly,  the  convergence  rate  increases  when  a  decreases.  This  is  consistent  with  the 
analytical  solution  shown  in  eqn.  (17)  for  k=l/a2,  where  the  convergence  rate  (f"ak/2)  is 
dependant  on  1/a  .  However,  when  a  is  too  small,  the  algorithm  becomes  unstable.  Figure  3-7 
shows  how  w  affects  the  convergence  rate  while  a  is  fixed  to  be  0.001.  Apparently,  the 
convergence  rate  increases  with  w .  When  w  is  too  large,  the  algorithm  again  becomes 
unstable. 

Figure  3-8  and  Figure  3-9  show  the  results  of  the  double  hump  model.  Clearly,  the 
extremum  seeking  algorithm  drives  the  cost  function  f  to  the  local  minimum. 

3.2  Vibration  Control  Testbed  Setup 

Figure  3-10  shows  a  detailed  sketch  of  the  whole  vibration  control  testbed  setup.  A  thin 
aluminum  cantilever  beam  with  one  piezoceramic  plate  bonded  to  each  side  is  fixed  on  a  block 
base  and  connected  to  the  ground.  The  two  piezoceramic  plates  are  used  to  excite  the  beam  by 
applying  electrical  field  across  their  thickness.  The  piezoceramic  plate  bonded  to  the  upper  side 
of  the  beam  is  called  the  "disturbance  piezoceramic"  because  it  is  used  to  apply  a  disturbance 
excitation  to  the  beam.  The  piezoceramic  plate  bonded  to  the  lower  side  of  the  beam  is  called 
the  "control  piezoceramic"  because  it  is  supplied  with  the  controller  output  signal  to  counteract 
the  disturbance  actuator.  The  beam  system  has  a  natural  frequency  of  about  97  Hz. 
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The  goal  of  the  disturbance  rejection  controller  is  to  mitigate  the  vibration  of  the 
aluminum  beam  generated  by  an  unknown  disturbance  signal.  The  controller  tries  to  generate  a 
signal  to  counteract  the  vibration  of  the  aluminum  beam  generated  by  the  “disturbance 
piezoceramic”.  The  performance  (or  the  residue)  signal  of  the  controller  is  measured  at  the 
center  of  the  tip  of  the  beam  by  a  laser-optical  displacement  sensor.  The  performance  signal  is 
filtered  by  a  high  pass  filter  with  fc=lHz  to  filter  out  the  dc  offset  of  the  displacement  sensor  and 

then  amplified  by  an  amplifier  with  a  gain  of  10. 

The  disturbance  and  control  signals  are  generated  by  our  dSPACE  (Model  DS1005)  DSP 
system  with  466MHz  PowerPC  CPU  and  amplified  by  two  separate  channels  of  an  amplifier  by 
a  same  gain  of  50.  The  types  and  conditions  of  the  signals  are  discussed  in  details  in  the  next 
section.  The  dSPACE  system  has  a  16-bit  A/D  and  a  16-bit  D/A  board.  The  computer  can 
acquire  data  using  Mlib/Mtrace  programs  in  MATLAB  through  the  dSPACE  system. 

The  whole  system  was  a  two-input/two-output  system.  One  input  was  the  control  signal 
and  the  other  input  was  the  “unmeasured”  disturbance  signal.  The  two  outputs  are  termed  a 
“reference  output”  and  a  “performance  output”.  For  this  validation  test,  the  reference  and 
performance  outputs  were  identical.  The  disturbance  rejection  algorithm  was  implemented  in  the 
Simulink  environment  and  compiled  and  downloaded  to  the  dSPACE  system.  The  disturbance 
signal  was  band-limited  white  noise  with  frequency  of  0-150  Hz. 

The  disturbance  rejection  algorithm  runs  in  one  of  the  following  two  modes:  1)  “ID,  then 
control”  (shown  in  Figure  3-11):  the  system  (control  model)  is  identified  by  the 
ARMARKOV/Toeplitz  system  ID  algorithm  and  the  identified  system  weight  matrix  Bzu  is 

transferred  to  the  ARMARKOV  control  algorithm;  then  the  controller  is  turned  on  and  the 
“control  signal”  is  switched  to  the  controller  output.  2)  “ID  and  control”  (shown  in  Figure  3-12): 
the  ID  and  control  are  turned  on  simultaneously.  The  input  (uID)  to  the  system  for  ID  can  be 

either  band-limited  white  noise  or  a  repetitive  linear  chirp  signal.  The  controller  uses  the 
identified  system  to  achieve  maximum  suppression  of  the  vibration  of  the  beam,  subject  to 
constraints  on  the  maximum  allowable  actuator  signal.  The  “ID  and  control”  mode  is  better 
when  the  system  is  a  time  variant  system  because  this  mode  updates  the  system  information 
during  every  iteration.  However,  the  “ID  and  control”  mode  adds  an  additional  signal  uID  to  the 

control  signal  all  the  time  and  this  certainly  affects  the  performance  of  the  disturbance  rejection 
controller.  The  tradeoff  between  the  adaptation  ability  and  effects  on  the  performance  should  be 
kept  in  mind. 

3.3  Results  of  the  Vibration  Control  Tests 
3.3.1  Computational  Tests 

For  real-time  control  applications,  the  turnaround  time  (defined  as  the  time  for  the 
program  to  execute  one  iteration)  is  required  to  be  less  than  the  sampling  time.  Complex 
algorithms  are  computationally  intensive  and  have  large  turnaround  time,  which  requires 
choosing  a  corresponding  larger  sampling  time  (or  a  smaller  sampling  frequency  fs ).  From  the 

Shannon  sampling  theorem,  the  sampling  frequency  must  be  larger  than  twice  the  highest 
frequency  of  interest  to  avoid  aliasing.  Thus,  algorithms  with  high  computational  complexity 
may  not  be  feasible  in  flow  control  applications.  The  tradeoff  between  choosing  a  large  fs  to 

satisfy  the  sampling  theorem  and  choosing  a  small  fs  to  allow  a  large  turnaround  time  must  be 
considered.  This  section  analyzes  the  effects  of  varying  the  parameters  of  the  ID  and  control 
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algorithms  on  the  computational  intensity.  This  serves  as  a  reference  for  choosing  the 
parameters  with  regard  to  the  computational  intensity. 

The  sampling  frequency  was  1024  Hz  for  the  computational  tests.  In  Figure  3-13  and 
Figure  3-14,  the  turnaround  time  of  the  system  ID  algorithm  by  varying  either  p,  n  or  p  is 
plotted,  while  the  other  two  parameters  are  fixed  at  unity.  It  is  shown  that  the  turnaround  time  is 
approximately  linearly  proportional  to  both  p  and  n ,  while  the  slope  for  n  is  approximately 
twice  of  that  for  p.  The  dependence  of  the  turnaround  time  on  p  is  approximately  quadratic. 
Clearly  the  averaging  window  number  p  has  much  more  significant  impact  on  the  computational 
intensity  than  the  other  two  parameters.  Figure  3-15  investigates  the  effects  of  varying  p  on  the 
computational  intensity  with  respect  to  n+p.  As  shown,  the  computational  intensity  is 
proportional  to  p.  From  these  results,  it  is  suggested  to  hold  p  to  be  a  small  number  and 
increase  p  to  improve  the  system  ID  performance. 

Figure  3-16  shows  the  effects  of  varying  nc  or  pc  on  the  computational  intensity  of  the 
ARMARKOV  disturbance  rejection  algorithm.  The  results  are  similar  to  those  of  the  system  ID 
algorithm.  The  turnaround  time  is  approximately  linearly  proportional  to  nc  and  pc  while  the 

slope  for  nc  is  approximately  twice  that  for  pc . 

3.3.2  System  Identification 

The  ARMARKOV/Toeplitz  system  ID  algorithm  requires  the  following  three  parameters: 
the  order  of  the  system  n ,  the  number  of  Markov  parameters  p ,  and  a  parameter  p  that 
determines  the  size  of  the  averaging  window.  The  SNR  is  also  a  parameter  that  can  affect  the 
performance  of  the  system  ID  algorithm.  When  the  number  of  Markov  parameters  is  unity,  the 
ARMARKOV  model  reduces  to  an  ARMA  model.  The  values  of  the  parameters  are  limited  by 
the  requirement  that  the  turnaround  time  must  be  less  the  sampling  time.  As  shown  in  the  last 
section,  p  has  the  smallest  effect  on  the  turnaround  time;  thus  in  this  section  the  performance  of 
the  system  ID  algorithm  with  varying  p  is  compared. 

The  offline  non-parametric  fit  of  the  frequency  response  of  the  beam  system  is  also 
implemented  as  a  comparison  and  shown  as  green  dot  lines  in  Figure  3-20  -  Figure  3-23.  The 
non-parametric  fit  uses  the  “invfreqz”  command  in  MATLAB  and  implements  a  second  order 
approximation.  The  “invfreqz”  command  returns  the  system  matrices  A  and  B  of  the  state 
space  representation.  The  zero-pole  map  of  the  non-parametric  fitted  system  is  shown  in  Figure 
3-17.  As  shown,  the  beam  system  is  a  low  damping  system  because  it  has  two  poles  that  are 
very  close  to  the  unit  circle.  The  controllability  matrix  is  [B  AB]  =  [l  1.655;0  l]  (for  a  second 

order  system),  which  has  full  rank  2.  This  means  that  the  system  is  controllable. 

The  sampling  frequency  was  1024  Hz.  The  input  signal  used  for  the  system  ID  was  a 
periodic  chirp  signal.  The  frequency  response  shown  in  Figure  3-20  -  Figure  3-23  was 
implemented  with  NFFT=1024,  no  overlap,  and  a  rectangular  window. 

Figure  3- 1 8  shows  a  very  good  match  between  the  measured  and  fitted  outputs  of  the 
system  with  the  system  ID  parameters  of  p=n=2,  p=10  and  SNR=20  dB.  Meanwhile,  as  shown 
in  Figure  3-19,  the  weight  tracks  of  the  system  ID  converge  at  about  0.5  seconds. 

Figure  3-20  to  Figure  3-23  show  the  comparison  between  the  measured  and  fitted 
frequency  response  with  the  system  ID  parameters  of  p=n=2 ,  SNR=20  dB  and  varying  p .  A 
significant  improvement  of  the  system  ID  is  obtained  when  p  changes  from  1  (ARMA  case)  to 
10.  Figure  3-24  compares  the  mean  square  error  (MSE)  verse  time  of  the  system  ID  with 
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varying  p.  Surprisingly,  it  is  found  that  the  case  with  p=10  has  the  best  performance.  This 
indicates  that  for  the  beam  system,  increasing  the  number  of  Markov  parameters  does  not 
necessarily  improve  the  system  ID  performance.  This  also  suggests  that  for  a  certain  system, 
there  may  exist  an  optimal  number  of  Markov  parameters.  It  is  also  shown  in  Figure  3-24  that 
for  the  case  with  p=40,  the  convergence  rate  is  slower  than  the  case  with  p=l ,  although  the  final 
MSE  is  better.  This  indicates  that  too  many  Markov  parameters  may  be  detrimental  to  the 
performance  of  the  ARMARKOV/Toeplitz  system  ID. 

Figure  3-25  compares  the  mean  square  error  (MSE)  verse  time  of  the  system  ID  with 
varying  SNR.  The  SNR  is  computed  by  using  the  formula:  SNR=101ogl0(o2/o2 ) ,  where  o2  is 

the  variance  of  the  control  signal  and  o2  is  the  variance  of  the  disturbance  signal.  It  is  clear  that 

the  system  ID  performs  better  with  a  higher  SNR.  This  suggests  that  when  the  system  ID  is 
conducted  with  unknown  disturbance,  it  is  better  to  apply  a  large  system  input  within  the 
maximum  allowable  range. 

3.3,3  Adaptive  Disturbance  Rejection 

The  ARMARKOV  disturbance  rejection  algorithm  requires  the  following  three 
parameters:  the  order  of  the  controller  nc ,  the  number  of  Markov  parameters  of  the  controller  pc 
and  the  adaptive  step  size  constant  y  that  controls  the  convergence  rate  of  the  controller. 

The  controller  uses  the  system  information  identified  by  the  ARMARKOV/Toeplitz 
system  ID  with  the  parameters  of  p=n=2,  p=10  and  SNR=20  dB  (shown  in  Figure  3-21).  All 
cases  of  the  controller  design  use  the  same  identified  system. 

Band-limited  white  noise  with  frequency  of  0  -  150  Hz  and  variance  of  0.09  is  used  to 
excite  the  disturbance  piezoceramics.  Figure  3-26  and  Figure  3-27  show  the  time  data  of  the 
performance  and  control  signals  with  the  parameters  of  nc=2 ,  pc=20  and  T  =  1 .  The  system  ID 

is  off  for  the  whole  period  and  the  controller  is  off  initially  and  turned  on  at  t  =  20  sec.  In  Figure 
3-28,  the  power  spectra  of  the  performance  signal  with  control  off  and  on  are  compared.  The 
power  spectra  were  calculated  by  using  the  time  data  of  20-second  duration  with  NFFT=1024, 
50%  overlap  and  a  hanning  window.  The  “control  on”  case  is  taken  after  the  controller  is  turned 

on  for  30  seconds.  The  performance  of  suppression  is  calculated  by  10glogl0  (a2ff  /o2n )  and  this 

case  gives  1 1.7  dB  suppression.  Interestingly,  in  Figure  3-28,  the  power  around  70  Hz  and  120 
Hz  of  the  “control  on”  case  is  higher  than  that  of  the  “control  off’  case.  This  is  generally  defined 
as  spillover  (Hong  and  Bernstein  1998).  Hong  and  Bernstein  (1998)  used  the  Bode  integral 
constraint  to  analyze  the  spillover  problem  and  concluded  that  the  spillover  is  inevitable  if  the 
reference  and  performance  signals  are  collocated  or  the  disturbance  and  control  actuators  are 
collocated.  For  this  vibration  control  test,  the  reference  and  performance  signals  are  collocated, 
thus  the  spillover  is  unavoidable. 

Figure  3-29,  Figure  3-30  and  Figure  3-32  show  the  performance  of  the  disturbance 
rejection  algorithm  with  varying  nc ,  pc  and  T ,  respectively,  while  other  parameters  are  held 
constant.  From 

Table  3-2  and  Table  3-3,  it  is  interesting  to  find  that  there  is  not  much  difference  of  the 
suppression  performance  for  varying  nc  and  pc .  However,  the  step  size  parameter  T  does  play 

a  significant  role.  Larger  T  gives  much  faster  convergence  and  better  performance.  However, 
if  T  is  too  large,  it  is  possible  for  the  controller  to  become  unstable.  This  tradeoff  should  be 
kept  in  mind  when  choosing  T  . 
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The  disturbance  rejection  controller  can  also  be  run  at  the  “ID  and  control”  mode.  In  this 
mode,  the  band-limited  white  noise  with  frequency  of  0  -  150  Hz  and  variance  of  0.09  is  used  to 
excite  the  disturbance  piezoceramics.  Meanwhile,  the  band-limited  white  noise  with  frequency 
of  0  -  150  Hz  and  variance  of  0.01  is  added  to  the  controller  output.  The  control  signal  is  shown 
in  Figure  3-33  and  the  performance  signal  is  shown  in  Figure  3-34.  Comparing  with  Figure 
3-27,  the  control  signal  of  the  “ID  and  control”  case  is  significantly  larger  at  the  beginning 
because  of  the  additive  signal  uID  and  the  evolution  of  the  controller  output  is  “buried”  under  it. 

Figure  3-35  compares  the  power  spectra  of  the  performance  signal  of  the  two  different 
modes.  It  is  surprising  that  the  “ID  and  control”  mode  results  in  lower  power  around  the  natural 
frequency  of  the  beam.  This  is  hard  to  see  in  Figure  3-34  because  the  “ID,  then  control”  mode 
seems  much  better.  However,  it  is  not  surprising  that  the  “ID  and  control”  mode  results  in  higher 
power  at  other  frequencies  than  the  natural  frequency  because  of  the  additive  signal  u1D.  In 
addition. 

Table  3-5  shows  that  the  “ID,  then  control”  mode  gives  better  suppression  performance 
of  the  overall  power. 

Unfortunately,  for  this  setup,  it  is  not  feasible  to  test  the  adaptability  of  the  two  modes. 
However,  this  will  be  done  in  the  wind  tunnel  experiments. 

As  a  summary,  the  computational  tests  are  conducted  first  to  determine  how  the 
parameters  affect  the  computational  complexity  of  the  system  ID  and  control  algorithms.  It  is 
shown  that  the  averaging  window  number  p  has  much  more  significant  impact  on  the 
computational  intensity  than  the  other  two  parameters  for  the  system  ID  algorithm.  The 
dependence  of  the  computational  complexity  vs.  n  is  approximately  twice  of  that  for  fi . 
Similarly  it  is  found  that  the  turnaround  time  of  the  control  algorithm  is  approximately  linearly 
proportional  to  nc  and  pc ,  while  the  slope  for  nc  is  approximately  twice  of  that  for  pc . 

The  ARMARKOV/Toeplitz  system  ID  algorithm  successfully  identifies  the  system 
(control  model)  and  results  in  very  good  frequency  response  approximations.  A  significant 
improvement  of  the  performance  of  the  ARMARKOV  system  ID  over  the  ARM  A  (when  p=l ) 
system  ID  is  found.  However,  too  many  Markov  parameters  of  the  ARMARKOV  system  ID 
may  be  detrimental  to  the  performance.  Higher  SNR  improves  the  performance,  thus  when  the 
system  ID  is  conducted  with  unknown  noise,  the  input  signal  should  be  chosen  as  large  as 
possible  within  the  maximum  allowable  level. 

The  order  of  the  controller  nc  and  the  number  of  Markov  parameters  pc  do  not  play 

significant  roles  on  the  performance  of  the  ARMARKOV  controller  for  this  vibration  control 
test.  However,  this  conclusion  may  vary  with  different  systems  and  remains  to  be  investigated. 
The  step  size  constant  T  significantly  affects  the  convergence  rate  of  the  controller.  T  should 
be  chosen  as  large  as  possible  before  it  makes  the  controller  unstable.  The  spillover  effect  is 
identified  in  this  vibration  control  test.  This  effect  is  unavoidable  because  the  “reference  signal” 
and  “performance  signal”  are  collocated  (Hong  and  Bernstein  1998). 

Table  3-1.  Parameters  for  the  simulations. 


Fs  (Hz) 

500 

Perturbation  amplitude 

a  =0.001,0.002,0.005 

Adaptation  gain 

k=l/a2 

Perturbation  frequency  ( Hz ) 

w  =  30,  40,  50 
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High  pass  filter  cutoff  frequency  ( Hz ) 

wh  =  10 

Low  pass  filter  cutoff  frequency  ( Hz ) 

w,=10 

Table  3-2.  Suppression  performance  of  the  disturbance  rejection  algorithm  with  pc=20,  Y  =  0.1 

and  varying  nc . 

pc  =20 ,  Y  =  0. 1  nc=l  nc=5 

nc  =  10 

Suppression  (dB)  8.9  8.6 

9.0 

Table  3-3.  Suppression  performance  of  the  disturbance  rejection  algorithm  with  nc-2,  Y  =  0. 1 


and  varying  juc . 

nc=2,Y  =  0.1  pc=l 

He  =20 

pc=40 

Suppression  (dB)  8.4 

8.7 

7.9 

Table  3-4.  Suppression  performance  of  the  disturbance  rejection  algorithm  with  nc=2,  pc=20 
and  varying  Y . 


n  =2 ,  li  =20  Y  =  0.01  Y  =  0.1  Y  =  1 

c  ’  r  c 


Suppression  (dB)  3.9  8.9  11.7 


Table  3-5.  Suppression  performance  of  the  disturbance  rejection  algorithm  at  different  modes 
with  nc  =2 ,  pc=20  and  Y  =  0.1. 


nc=2,pc=20  and  Y  =  0.1 

ID,  then  control 

ID  and  control 

Suppression  (dB) 

7.1525 

4.9994 
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Cost  function 


20: - ; - r- 


Figure  3-1.  One-dimensional  example  of  the  downhill  simplex  algorithm. 


Figure  3-2.  Two-dimensional  example  of  the  downhill  simplex  algorithm. 
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Figure  3-3.  Simulation  block  diagram  for  extremum  seeking  control. 
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Figure  3-4.  Single  global  maximum  test  model:  f=f*-(0-0”)  ,  where  f *=10  and  0*=5. 
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Figure  3-5.  Double  hump  model  with  one  local  maximum  and  one  global  maximum.  The 


^  Jf  ,  .  f  f  (0)=-1.2e'1208-1.6e,o07  +  1.4e'706-2.3e'505 

function  is  fitted  by  a  polynomial: 

+  1.5e'304-4.4e‘203+3.8e‘,02+2.30+5.4 


o 

a=0.002 

;  '  •  '  '  i 

a=0.005 

.51 _ i _ ; _ _ _ j _ i - - - i — .  r  ,.i  — ! 

0  10  2030405060  70  80  90  100 

t  (sec) 


a=0.005 

^0  10  20  30  40  50  GO  70  80  90  100 

t  (sec) 

Figure  3-6.  6  converges  to  the  optimal  input  0*=5  and  f  converges  to  the  global  maximum 
f*  =  10 at  various  convergence  rates  due  to  various  values  of  a  while  w=50. 
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Figure  3-7. 


Figure  3-8. 
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0  converges  to  the  optimal  input  0*=5  and  f  converges  to  the  global  maximum 
f  =10  at  various  convergence  rates  due  to  various  values  of  w  while  a=0.001 . 


-2 

0  50  100  150  200  250  300  350  400 

t  (sec) 


0  converges  to  the  local  optimal  input  (see  Figure  3-5)  0*  =  14  .  (a=0.001,w=50). 
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Figure  3-9.  f  converges  to  the  local  maximum  (see  Figure  3-5)  f=40.  (a=0.001,w=50). 


48 


Performance  Signal 


PC 

ML  IB/M  1  KAC  L 

!  A 

dSPACt 

DS  1005 

C  h2  9  CM 


DS200I 
A/D  Board 


DS2102 
D/A  Board 


\  th  2 


High  Pass  Filter 
and  Amplifier 


Aluminum  Shim 


Control  _ 
Pic/oceramic 


OptoNCDT  2000 
Laser-Optical 
Displacement  Sensor 


Disturbance 

Piezoceramic 


Figure  3-10.  Vibration  Control  Testbed. 


white  noise  or 
chirp 

Figure  3-11.  Block  diagram  of  vibration  control  with  “ID,  then  control”. 
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while  noise  or 
chirp 

Figure  3-12.  Block  diagram  of  vibration  control  with  “ID  and  control”. 
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Figure  3-13.  Effects  of  varying  n,  n  or  p 
ARMARKOV/Toeplitz  system  ID. 
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Figure  3-14.  Effects  of  varying  p  on  the  computational  intensity  of  the  ARMARKOV/Toeplitz 
system  ID. 
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Figure  3-15.  Effects  of  varying  p  on  the  growth  rate  of  the  computational  intensity  of  the 
ARMARKOV/Toeplitz  system  ID  with  respect  to  n+p . 
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Figure  3-16.  Effects  of  varying  nc  or  pc  on  the  computational  intensity  of  the  ARMARKOV 
disturbance  rejection. 


1  j 

J 

0.8  i 

iT\ 

i 

•j 

0.6 ! 

0.4' 

j 

i 

tr 

05 

CL 

CVJ 

o 

[ 

03 

C 

0 

i 

j 

CT) 

05 

E 

-0.2  ^ 

-0.4 

1 

j 

-0.6; 

-[ 

-0.8 ' 

i 

i 

i 

-i  \ 

_ l _ j 

-i 

-0.5  0  0.5 

1 

Real  Part 

Figure  3-17.  Zero-pole  map  of  the  non-parametric  fit  of  the  frequency  response. 
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Figure  3-18.  Measured  output  and  fitted  output  by  the  ARMARKOV/Toeplitz  system  ID  with 
p=n=2 ,  p=10  and  SNR=20  dB. 


0.8 

n 


0.7 


|  06  j|V  j j 

0.5 1  "j 


2  3  4  5  6  7 

t  (sec) 


10 


Figure  3-19.  Weight  tracks  of  the  ARMARKOV/Toeplitz  system  ID  with  p=n=2,  p=10  and 
SNR=20  dB. 
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Figure  3-20.  Comparison  of  measured  frequency  response,  non-parametric  fit  and  the 
ARMARKOV/Toeplitz  system  ID  with  p=n=2,  p=l  and  SNR=20  dB. 
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Figure  3-21.  Comparison  of  measured  frequency  response,  non-parametric  fit  and  the 
ARMARKOV/Toeplitz  system  ID  with  p=n=2,  p=10  and  SNR=20  dB. 
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Figure  3-22.  Comparison  of  measured  frequency  response,  non-parametric  fit  and  the 
ARMARKOV/Toeplitz  system  ID  with  p=n=2,  p=20  and  SNR=20  dB. 
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Figure  3-23.  Comparison  of  measured  frequency  response,  non-parametric  fit  and  the 
ARMARKOV/Toeplitz  system  ID  with  p=n=2,  p=30  and  SNR=20  dB. 
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Figure  3-24.  Comparison  of  MSE  of  the  ARMARKOV/Toeplitz  system  ID  with  p=n=2, 
SNR=20  dB  with  different  p . 
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Figure  3-25.  Comparison  of  MSE  of  the  ARMARKOV/Toeplitz  system  ID  with  p=n=2, 
p=20with  different  SNR. 
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Figure  3-26.  Performance  signal  of  the  ARMARKOV  disturbance  rejection  to  band-limited 
white  noise  (0-150  Hz)  with  nc=2,  p.c  =20  and  Y  =  l. 
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Figure  3-27.  Control  signal  of  the  ARMARKOV  disturbance  rejection  with  nc=2,  nc=20  and 
T  =  1 . 
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Figure  3-28.  Power  spectra  of  the  performance  signals  of  the  ARMARKOV  disturbance 
rejection  with  nc=2,  pc=20  and  Y  =  l. 


20 - r  - 

Control  off 

Control  on  n  =1  ! 

c 

q  Control  on  nc= 5  [ 

Control  on  n  =10 
c 


Frequency  (Hz) 

Figure  3-29.  Comparison  of  power  spectra  of  the  performance  signals  of  the  ARMARKOV 
disturbance  rejection  with  T  =  0.1 ,  pc=20  and  different  nc . 
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Figure  3-30.  Comparison  of  power  spectra  of  the  performance  signals  of  the  ARMARKOV 
disturbance  rejection  with  nc=2 ,  T  =  0.1  and  different  pc . 
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Figure  3-31.  Comparison  of  convergence  of  the  ARMARKOV  disturbance  rejection  with  nc=2  , 
pc=20  and  different  T  . 
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Figure  3-32.  Comparison  of  power  spectra  of  the  performance  signals  of  the  ARMARKOV 
disturbance  rejection  with  nc=2,  pc=20  and  different  T. 


i.f  '  il‘ 


I  .  •  |  ,  y  '  t 

-0.3 :  - 

°40  5  10  15  20  25  30  35  40 

t  (sec) 

Figure  3-33.  Control  signal  of  the  ARMARKOV  disturbance  rejection  at  “ID  and  control”  mode 
with  nc=2,  pc=20  and  Y  =  0.1. 
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Figure  3-34.  Comparison  of  convergence  of  the  ARMARKOV  disturbance  rejection  at  different 
modes  with  nc=2,  pc=20  and  Y  =  0.1. 
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Figure  3-35.  Comparison  of  power  spectra  of  the  performance  signals  of  the  ARMARKOV 
disturbance  rejection  at  different  modes  with  nc=2 ,  ptc  =20  and  Y  =  0.1 . 
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4  Experimental  Setup  and  Data  Analysis  Method 

The  separation  control  experiments  are  conducted  in  an  open-return  low-speed  wind 
tunnel  with  a  30.48  cm  (1  ft)  by  30.48  cm  test  section.  The  wind  tunnel  has  two  anti-turbulence 
screens,  an  aluminum  honeycomb  and  a  9:1  contraction  ratio.  The  airspeed  is  controlled  by  the 
variable  frequency  of  the  motor  fan. 

A  two-dimensional  NACA  0025  airfoil  that  is  equipped  with  synthetic  jet  actuators, 
Kulite  dynamic  pressure  transducers  and  a  lift/drag  balance  is  used  as  the  test  model.  A  Particle 
Image  Velocimetry  (PI V)  system  is  used  for  flow  visualization  and  quantitative  flow  field 
measurements.  A  Dantec  CTA  hot  wire  system  is  used  to  measure  instantaneous  velocity. 

This  chapter  describes  each  part  of  the  experimental  setup  in  detail.  A  brief  description 
of  the  Higher  Order  Statistical  Analysis  (HOSA)  is  also  presented  in  this  chapter  because  it  may 
be  used  for  nonlinear  flow  instability  analysis. 

4.1  NACA  0025  Airfoil  Model 

A  two-dimensional  NACA  0025  airfoil  with  chord  length  of  15.24  cm  (6  in.)  is  built  as  a 
test  bed  for  flow  separation  control  (Figure  4-1).  The  span  of  the  airfoil  model  is  29.21  cm  (1 1.5 
in.),  which  allows  for  a  slight  gap  on  either  end  to  accommodate  a  sidewall-mounted  strain- 
gauge  sting  balance.  The  boundary  layer  is  tripped  at  the  leading  edge  region  using  No.  60  sand 
grit.  Two  pairs  of  synthetic  jets  are  embedded  in  the  airfoil  at  approximately  3%  chord  and  30% 
chord,  respectively.  Six  ports  near  the  rear  of  the  airfoil  at  the  mid-span  location  are  available 
for  dynamic  pressure  transducers.  The  six  ports  are  located  at  approximately  44.0%,  52.5%, 
61.0%,  69.5%,  77.9%  and  86.4%  chord.  A  pre-amplifier  PCB  board  for  the  dynamic  pressure 
transducers  can  be  also  installed  in  the  airfoil.  The  detailed  side  view  of  the  airfoil  is  shown  in 
Figure  4-1. 

4.2  Synthetic  Jet  Actuators 

The  airfoil  is  fitted  with  two  pairs  of  synthetic  jet  arrays  (each  with  0.5  mm  wide  slots 
separated  by  2.4  mm),  which  are  located  in  the  central  l/3rd  spanwise  region  of  the  airfoil  (see 
Figure  4-1).  The  first  pair  is  located  near  the  leading  edge  of  the  airfoil,  at  approximately  3% 
chord,  while  the  second  is  placed  near  the  point  of  maximum  thickness  at  about  30%  chord.  The 
first  array  is  fixed,  while  the  second  array  can  be  translated  between  25%  chord  and  37%  chord. 
The  detailed  design  procedures  of  the  synthetic  jet  actuators  that  are  used  in  this  research  can  be 
found  in  Gallas  et  al.  (2003)  and  Gallas  (2005).  The  primary  goal  of  the  design  is  to  maximize 
the  magnitude  of  the  volume  flow  rate  through  the  orifice  per  applied  voltage  (i.e.  |Qoul/Vac|, 

where  Qoul  is  the  volume  flow  rate  and  is  the  applied  ac  voltage)  over  a  frequency  range  of 

O(kHz),  while  the  size  of  the  synthetic  jet  actuators  are  limited  by  the  geometry  of  the  airfoil. 

As  mentioned,  the  frequency  response  of  the  synthetic  jet  actuators  is  another  important  design 
criterion.  The  frequency  response  of  the  actuators  must  be  chosen  appropriately  to  effectively 
control  (via  amplitude  and  burst  modulation  techniques)  the  flow  separation  over  a  range  of 
frequencies,  ranging  from  the  low  frequency  shedding  in  the  wake  to  the  high  frequency  shear 
layer  instability.  The  side  and  top  views  of  the  synthetic  jet  actuators  are  shown  in  Figure  4-3. 
The  cavity  is  151  mm  long,  28  mm  high  and  2  mm  wide.  Five  piezoceramic  disks  are  attached 
to  one  side  of  the  cavity.  They  are  driven  in  phase  using  a  single  amplified  drive  signal  to 
achieve  maximum  flow  rate.  A  thin  slot  (0.5  mm  wide  by  101.6  mm  long)  at  the  top  of  the 
cavity  permits  oscillatory  fluid  flow.  Two  closely  spaced  synthetic  jets  are  obtained  by 
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introducing  a  rigid  wall  to  separate  them,  as  shown  in  Figure  4-3.  The  two  synthetic  jet  actuators 
are  nominally  identical.  See  the  detailed  characterization  of  the  actuators  in  0. 

4.3  Experimental  Methods 

4.3.1  Flow  Visualization 

Flow  field  velocity  data  over  the  surface  of  the  airfoil  and  in  the  wake  are  acquired  using 
Particle  Image  Velocimetry  (PIV).  The  PIV  system  consists  of  a  pair  of  New  Wave  Minilase  15 
Hz,  50  mJ  per  pulse,  Nd:YAG  lasers  with  appropriate  light  sheet  optics.  The  width  of  the  light 
sheet  is  approximately  1-2  mm  at  the  plane  of  measurement.  A  TSI  model  630157  Powerview 
Plus  2MP  10-bit  CCD  camera  is  used  to  acquire  images.  This  camera  contains  1600  x  1200,  7.4 
pm  square  pixels.  A  series  of  Nikon  lenses  (60  mm,  75-240  mm,  200  mm)  are  available.  The 
flow  is  seeded  with  water-based  fog  fluid  by  a  LeMaitre  G150  seeder  and  the  seeding  density  is 
adjusted  to  insure  uniform  seeding  density. 

The  laser  pulse  generator  and  the  camera  are  synchronized  by  a  TSI  Model  610032 
Synchronizer  which  is  configured  to  acquire  a  pair  of  images  using  TSI  INSIGHT  Software 
version  6.1.1.  The  computation  of  the  velocity  field  begins  by  dividing  the  image  into  a  grid  of 
interrogation  windows  overlapped  in  space  by  50%.  These  windows  typically  range  from  32  x 
32  pixels  to  64  x  64  pixels.  The  velocity  is  determined  by  the  known  distance  that  a  particle  is 
displaced  during  the  known  time  dT.  The  INSIGHT  Software  utilizes  an  FFT  cross-correlation 
process  in  conjunction  with  a  Gaussian  peak  search  algorithm  to  calculate  the  average  velocity  of 
the  particles  in  the  interrogation  widow.  A  number  of  validation  schemes  are  available  in  the 
software,  such  as  range  outlier  rejection  and  median  filtering. 

4.3.2  Lift/Drag  Balance 

A  strain-gauge  balance  is  designed  to  measure  lift  and  drag  forces  of  the  airfoil  test  bed. 
The  detailed  design  procedure  can  be  found  in  Griffin  (2003).  Two  pairs  of  strain  gauges  are 
attached  to  the  cantilever  that  supports  the  airfoil  to  measure  the  normal  and  axial  forces  on  the 
airfoil,  respectively  (Figure  4-4).  The  layout  of  the  strain  gauges  and  Wheatstone  bridge 
configuration  are  shown  in  Figure  4-5  and  Figure  4-6.  The  output  of  the  Wheatstone  can  be 
calculated  by  the  following  equation 


VAR 

R 


(81) 


From  the  above  equation,  we  can  see  that  the  output  is  linearly  dependant  on  the  change 
of  resistance  AR  . 

The  output  of  the  strain  gauges  is  measured  by  a  high-resolution  HP34970A  DAQ  system 
and  is  averaged  over  2  power  line  cycles  to  eliminate  60  Hz  noise.  The  lift  and  drag  are 
calculated  from  the  normal  and  axial  forces  together  with  the  angle  of  attack  (Figure  4-4)  via  the 
following  equations: 

L=Ncos(AOA)-Asin(AOA)  (82) 

D=Nsin(AOA)+Acos(AOA)  (83) 


where  L  and  D  stand  for  lift  and  drag,  respectively,  while  N  and  A  stand  for  normal  and  axial 
force,  respectively. 

Before  the  balance  can  be  used  for  the  wind  tunnel  experiments,  it  is  calibrated  by  adding 
known  weights  on  the  balance  and  measuring  the  output  from  the  normal  and  axial  strain  gauges. 
Figure  4-7  and  Figure  4-8  show  typical  normal  and  axial  force  calibrations  vs.  balance  output. 
Very  good  linear  relationships  between  the  balance  output  and  the  forces  on  the  balance  are 
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achieved.  The  coefficients  of  the  linear  equations  are  used  to  back  out  the  forces  on  the  airfoil 
from  the  voltage  output  of  the  strain  gauges. 

The  balance  is  also  validated  by  comparing  with  the  lift  and  pressure  drag  coefficients 
measured  by  integrating  the  static  pressure  around  the  airfoil.  The  pressure  taps  are  located  at 
the  center  span  of  the  airfoil  and  the  static  pressure  is  measured  via  a  Heise  static  pressure  gauge. 
Figure  4-9  shows  the  static  pressure  distributions  on  the  airfoil  surface  at  different  AOAs  when 
Re  =  150, 000.  From  Figure  4-9,  it  can  been  identified  that  the  flow  is  separated  at  AOA=15° 

and  20° .  The  suction  zones  on  the  upper  surface  shrink  dramatically.  This  is  generally  referred 
to  as  pressure  loss  due  to  flow  separation  and  is  primarily  responsible  for  deteriorating  lift  to 
drag  ratio. 

The  lift  and  drag  coefficients  are  calculated  by  integrating  the  static  pressure  around  the 
airfoil  surface,  assuming  surface  friction  is  negligible  compared  with  pressure  forces  and  the 
flow  is  two-dimensional.  Figure  4-10  and  Figure  4-11  show  the  comparison  of  the  lift  and  drag 
coefficients  calculated  by  the  two  different  methods  at  Re  =  100, 000  and  Re  =  150, 000.  The 

2o 

uncertainty  was  calculated  at  %95  confidence  interval  (i.e.  uncertainty  =  —j= ,  where  o  is 

V  N 

standard  deviation  and  N  is  number  of  measurements).  As  shown,  they  agree  reasonably  well 
considering  measurement  uncertainties.  This  validates  that  the  balance  works  as  desired.  The 
main  reason  for  the  differences  is  the  three-dimensional  effect  as  the  pressure  taps  only  measure 
at  the  center  span. 

4.3.3  Dynamic  Pressure  Transducers 

To  measure  the  pressure  fluctuation  on  the  airfoil  surface,  it  is  required  that  the  pressure 
sensors  must  be  compact  so  that  they  can  be  installed  within  the  limited  space  in  the  airfoil.  It  is 
also  desired  that  they  have  large  enough  bandwidth  to  capture  the  characteristics  of  the 
oscillations  of  the  flow  above  the  airfoil,  and  their  response  is  linear  with  respect  to  the  pressure 
load  within  the  range  of  interest.  For  these  reasons,  a  number  of  commercially  available  MEMS 
Kulite  LQ 125-5 A  dynamic  pressure  transducers  (Figure  4-12)  are  used  to  obtain  dynamic 
pressure  response  on  the  upper  surface  of  the  airfoil.  The  transducers  can  be  flush  mounted  in 
the  six  available  locations  on  the  upper  surface.  A  pre-amplifier/filter  board  for  the  transducers 
is  designed  to  eliminate  dc  response  (fcutoff  =  1.5  Hz )  and  amplify  the  outputs  by  a  gain  of  100. 

The  pre-amplifier/filter  board  can  be  installed  inside  the  airfoil  so  that  the  airfoil  acts  like  an 
electronic  enclosure.  Before  the  transducers  can  be  used  in  the  experiments,  they  are 
dynamically  calibrated  in  a  2.54  cm  (1  in.)  by  2.54  cm  plane  wave  tube  (PWT).  A  speaker  was 
used  as  a  source,  and  a  Briiel  &  Kjaer  (Model  4318)  microphone  was  used  as  a  reference 
transducer.  Figure  4-14  shows  the  linear  response  of  a  typical  Kulite  sensor  that  is  obtained  by 
fixing  the  frequency  and  increasing  the  input  amplitude  of  the  speaker.  The  frequency  response 
is  measured  using  a  periodic  chirp  signal  (Figure  4-14).  As  shown,  the  frequency  response  does 
not  vary  up  to  approximately  3000  Hz,  which  is  sufficient  for  this  research. 

4.3.4  Hot  Wire  Anemometry 

A  Dantec  constant-temperature  hot  wire  anemometry  system  (CTA  module  90C10)  is 
used  to  measure  time-resolved  velocity  in  the  unseparated  flow  above  the  airfoil.  The  CTA 
system  includes  A/D  converter  and  all  the  signal  conditioners  needed.  Before  the  measurements, 
a  static  calibration  is  performed  by  the  calibration  module  and  the  flow  unit  (90H01  and  90H02). 
A  typical  calibration  curve  is  shown  in  Figure  4-15.  Since  the  output  of  the  hot  wire  system 
usually  drifts  due  to  temperature  changes,  connections,  etc,  the  calibration  should  be  done  before 
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each  measurement.  Two  algorithms  are  commonly  used  for  curve  fitting.  One  is  a  polynomial 

that  is  used  here,  and  the  other  is  King’s  law  (power  law):  U=((E2-A)/b)  ,  where  E  is  the 

voltage  output  of  the  hot  wire  and  U  is  the  flow  velocity  (Jprgensen  1996).  The  difference 
between  the  temperatures  at  calibration  and  measurements  should  also  be  corrected  by  means  of 

EC0„=E,^-y  j  ,  where  Tw  is  the  wire  temperature,  T0  is  the  temperature  at  calibration. 


£j  is  the  raw  wire  voltage,  T,  is  the  temperature  during  measurement  and  Econ.  is  the  corrected 

voltage  (Jprgensen  1996).  During  experiments,  the  hot  wire  probe  (55P1 1)  is  mounted  on  a  2- 
dimensional  Velmex  traversing  system  which  has  spatial  resolution  of  about  1 ,6pm/step  in  both 
directions. 

4.4  Control  System  Hardware  and  Software 

The  control  systems  for  the  separation  control  experiments  are  implemented  by  a 
dSPACE  (Model  DS1005)  DSP  system  with  a  466MHz  PowerPC  CPU.  The  dSPACE  system 
has  a  5-channel  16-bit  A/D  board  (DS2001)  and  a  6-channel  16-bit  D/A  board  (DS2102)  as  the 
data  acquisition  equipments.  The  range  of  the  data  acquisition  boards  can  only  be  -10  to  +10  V, 
0  to  10  V  or  -5  to  +5V.  The  control  algorithms  are  first  programmed  in  Matlab/Simulink  and  C 
programs  (c-mex  sfunction)  and  then  compiled  and  downloaded  to  the  dSPACE  system.  The 
compiled  programs  together  with  the  data  acquisition  boards  are  able  to  run  the  experiments  in 
real  time.  The  computer  is  also  able  to  acquire  data  into  Matlab’s  workspace  through  the 
dSPACE  system  via  the  m-lib  programs  provided  by  the  dSPACE. 

4.5  Higher  Order  Statistical  Analysis  (HOSA) 

Higher  order  spectral  analysis  is  used  to  uncover  the  nonlinear  interactions  in  signals  or 
to  identify  nonlinear  systems  (Nikias  and  Mendel  1993).  As  discussed  in  Error!  Reference 
source  not  found.,  there  are  three  characteristic  frequencies  and  nonlinear  interactions  between 
them  are  inherent  in  separated  flow.  Unfortunately,  the  power  spectrum  alone  is  incapable  of 
providing  any  conclusive  proof  of  the  nonlinear  interactions.  The  power  spectrum  only  provides 
proof  of  presence  of  power  at  certain  frequencies.  On  the  other  hand,  higher-order  spectral 
method  can  quantify  quadratic  coupling  between  frequency  pairs.  For  example,  it  can  provide 
the  information  that  the  generation  of  power  at  a  certain  frequency  is  the  result  of  quadratic 
coupling  of  other  frequencies.  The  auto-bispectrum  uses  third  order  cumulants  and  is  defined  as 


(84) 


and  the  auto-bicoherence  is  defined  as 


B. 


(85) 


where  X(f)  denotes  the  Fourier  transform  of  x(t),  *  denotes  the  complex  conjugate  and 
Pxx  (f )  denotes  the  auto-spectrum  of  x(t). 

The  auto-bicoherence  is  bounded  by  zero  and  unity.  Disturbances  with  frequencies  f, ,  fj 
and  f +fj  are  quadratically  coupled  if  b2  (fi,fj)  =  l ,  not  quadratically  coupled  if  b:  (f,,fj)=0  and 
partially  coupled  if  0<b2  (fi,fj)<l . 
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Just  as  the  auto-spectrum  has  the  cross-spectrum  as  its  counterpart  for  signals  x(t)  and 
y(t),  the  auto-bi spectrum  has  the  cross-bispectrum  as  its  counterpart,  which  is  defined  as: 

B..,(f,.fj)  =  ltalf[x(fi)X(fj)Y-(f1+fJ)]  (86) 

Similarly,  the  cross-coherence  is  obtained  by  normalizing  the  cross-bispectrum  and 
defined  as  follows: 

|2 


bUCfih 


B„,  (f|.fj)| 


(87) 


Figure  4-1.  NACA  0025  airfoil  model  with  actuators  and  pressure  transducers  installed. 
(Adapted  from  Holman  et  al.  2003) 
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Figure  4-2.  Schematic  of  a  synthetic  jet  actuator. 
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Figure  4-3.  Synthetic  jet  array.  (Adapted  from  Holman  et  al.  2003) 
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Figure  4-4.  Forces  on  NACA0025  airfoil 
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Normal  force 


Figure  4-5.  Closer  view  of  the  strain  gauges. 


Figure  4-6.  Wheatstone  bridge  configuration  of  the  balance. 
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Figure  4-7.  Normal  force  vs.  balance  output 
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Figur  Axial  force  vs.  balance  output 
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Figure  4-9.  Static  pressure  distributions  on  the  airfoil  surface  at  different  AOA  at  Re  =  150,000 . 
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Figure  4-10.  Comparison  of  lift  and  drag  coefficients  measured  by  the  static  pressure  and  the 
balance  at  Re  =  100, 000. 
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Figure  4-11.  Comparison  of  lift  and  drag  coefficients  measured  by  the  static  pressure  and  the 
balance  at  Re  =  150, 000. 
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Figure  4-12.  Picture  of  a  Kulite  transducer. 
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Figure  4-13.  Linear  response  (at  500Hz)  of  a  typical  Kulite  transducer. 
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o  Frequency  response  of  a  typical  Kulite  transducer. 
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5  Results  and  Discussion 


The  experimental  results  of  the  dynamic  feedback  control  and  nonlinear  control 
approaches  are  presented  in  two  parts. 

5.1  Dynamic  Feedback  Control 

5.1.1  Experimental  Configuration 

Figure  5-1  shows  the  complete  experimental  configuration.  The  system  ID  and  control 
algorithms  run  on  the  dSPACE  controller  system  in  real-time  (~4  KHz).  The  controller 
generates  control  signal  that  is  amplified  by  an  amplifier.  The  pressure  fluctuation  signals 
measured  by  Kulite  sensors  are  amplified  and  filtered  before  sending  to  the  dSPACE  controller. 

5.1.2  System  Identification 

5.1.2.1  Coherent  flow  structures 

Unlike  POD-based  approaches  (Holmes  et  al.  1998;  Tadmor  et  al.  2007;  Ausseur  et  al. 
2007),  the  unsteady  surface  pressure  signals  are  used  exclusive  of  the  velocity  field  for  feedback 
in  this  research.  Although  system  ID  and  POD-based  methods  of  modeling  the  flow  are 
different,  they  all  attempt  to  capture  the  signature  of  the  separated  flow  -  the  coherent  flow 
structures.  In  this  section,  we  show  that  the  surface  pressure  signals  indeed  represent  the 
footprint  of  the  coherent  flow  structures. 

The  following  experiment  is  devised  to  show  this.  A  continuous  pulse  train  (repetition 
rate  at  1  Hz  and  amplitude  at  50  V)  is  fed  to  the  actuator  Al  at  Re=  120,000.  The  pressure 
signals  are  phase  averaged  relative  to  the  pulse  signal.  Four  thousand  averages  are  taken  to 
obtain  the  statistically  converged  pressure  signal  profiles.  As  shown  in  Figure  5-2,  a  vortex 
(produced  by  the  actuator  pulse)  propagates  downstream  indicated  by  the  surface  pressure 
fluctuations.  The  vortex  reaches  the  sensor  SI  first  and  then  S2  ...  S6.  After  the  vortex  passes 
by,  the  averaged  surface  pressure  fluctuating  approaches  zero.  This  is  because  the  random 
pressure  fluctuations  that  are  not  correlated  with  the  pulse  input  from  the  actuator  possess 
random  phase  and  are  averaged  out.  The  convective  velocity  is  much  slower  than  the  free  stream 
velocity,  with  a  nominal  lifetime  of  more  than  5  airfoil  chord  lengths.  With  the  link  between  the 
flow  structure  and  the  surface  pressure  clearly  established,  the  goal  of  the  system  ID  approach  is 
to  correlate  the  actuator  input  and  the  corresponding  surface  pressure  fluctuations  and  utilize  the 
relationship  to  model  the  coherent  flow  structure  with  linear  dynamical  equations. 

5. 1.2.2  Linear  prediction 

The  control  approach  is  based  on  the  assumption  that  the  coherent  flow  structures  may  be 
modeled  by  linear  dynamical  equations.  This  section  demonstrates  that  the  linear  model  is 
capable  of  predicting  the  downstream  evolution  of  the  flow  dynamics  (measured  by  the  pressure 
sensors)  subject  to  the  actuation  upstream  (provided  by  the  ZNMF  actuators).  In  fact,  previous 
studies  in  turbulent  boundary  layer  (Rathnasingham  and  Breuer  2003)  and  cavity  flows 
(Cattafesta  et  al.  1999)  have  suggested  that  linear  approximations  can  reasonably  predict 
inherently  nonlinear  flow  structures. 

The  system  ID  algorithm  is  applied  first  to  demonstrate  this.  The  computational 
requirements  are  demanding  for  both  the  system  ID  and  control  algorithms.  When  implementing 
the  algorithms,  the  algorithm  parameters  were  chosen  due  to  hardware  limitations  and  were  not 
optimized.  The  sampling  frequency  is  chosen  to  be  4096  Hz.  In  Figure  5-3,  z  denotes  the 
pressure  signal  measured  by  the  #6  pressure  sensor  shown  in  Figure  5-1  subject  to  a  band-limited 
random  input  provided  to  actuator  Al.  Using  the  input  and  z,  one  can  fit  a  model  to  represent 
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the  flow  structures  using  the  approach  described  in  earlier  sections.  One  should  be  aware  that 
this  model  actually  includes  actuator  and  sensor  dynamics  and  other  hardware  in  the  loop,  e.g. 
amplifiers  and  filters.  Here,  Error!  Objects  cannot  be  created  from  editing  field  codes,  in 
Figure  5-3  is  the  estimated  Error!  Objects  cannot  be  created  from  editing  field  codes,  using 
the  model  mentioned  above,  and  the  data  show  a  reasonable  match  over  a  range  of  time  or 
frequency  scales.  The  errors  are  primarily  due  to  the  random  turbulent  structures  that  are 
uncorrelated  with  the  actuator  input,  and  they  cannot  be  modeled  by  the  system  ID  algorithm. 

An  indication  of  the  convergence  of  the  system  ED  to  the  actual  flow  model  is  the 

expected  value  of  the  Mean  Squared  Error  (MSE),  which  is  defined  as  MSE=E^[z(k)-z(k)]‘  J. 

As  shown  in  Figure  5-4,  at  the  beginning  the  MSE  has  a  relatively  large  value.  This  is  because 
the  model  parameters  are  initialized  to  zero.  Then  the  model  parameters  are  trained  by  the 
ARMARKOV  system  ID  algorithm  targeting  the  objective  of  minimizing  the  error. 

5. 1.2.3  Frequency  response  and  the  performance  of  system  ID 

To  further  evaluate  the  performance  of  system  ID  approach,  we  compare  the  frequency 
response  of  the  flow  system  determined  using  conventional  FFT  methods  for  single-input/single¬ 
output  systems  described  in  Chapter  6  of  Bendat  and  Piersol  (2000)  with  that  determined  using 
the  converged  ARMARKOV  system  ID  model  parameters.  When  computing  the  frequency 
response,  the  parameters  are  fs=4096  Hz,  NFFT=1024,  75%  overlap,  a  Hanning  window  and 
320  effective  averages. 

Figure  5-5  shows  that  the  system  ID  does  not  perfectly  match  the  frequency  response,  but 
it  does  capture  the  essential  characteristics  over  a  broad  frequency  range.  It  is  also  clear  that  the 
coherence  between  the  input  and  output  is  close  to  zero  at  frequencies  less  than  600  Hz,  which  is 
a  characteristic  of  the  present  piezoelectric  zero-net  mass  flux  actuators,  which  posses  a 
resonance  near  1200  Hz  (Holman  et  al.  2003).  The  low  coherence  renders  the  FFT-based 
frequency  response  estimate  uncertain  and  highlights  the  difficulty  of  designing  a  control  system 
using  classical  frequency  domain  approaches. 

5. 1.2.4  Acoustic  contamination 

In  all  of  the  discussions  above,  we  have  ignored  an  important  potential  issue  related  to 
acoustic  contamination.  It  is  well  known  that  zero-net  mass  flux  actuators  can  produce 
significant  sound.  In  the  present  control  problem,  the  pressure  sensors  are  intended  to  capture 
the  hydrodynamics  of  the  coherent  flow  structures.  However,  as  demonstrated  by  Figure  5-6,  the 
pressure  sensors  do  not  discriminate  between  acoustic  and  hydrodynamic  pressure  fluctuations. 
Since  the  pressure  measurements  contain  both  components,  the  disturbance  rejection  control 
algorithm  will  try  to  suppress  the  acoustic  power  as  well  as  the  hydrodynamic  power,  possibly 
resulting  in  an  undesirable  reduction  in  the  actuator  amplitude.  Furthermore,  from  a  control 
standpoint,  the  acoustic  and  hydrodynamic  paths  have  significantly  different  propagation  speeds, 
leading  to  significant  phase  lag  differences  and  an  unstable  controller. 

One  way  to  address  this  problem  is  to  estimate  a  frequency-wavenumber  spectrum  using 
Fourier-based  methods,  but  such  a  method  is  not  amenable  to  a  real-time  control  system.  A 
second  approach,  adopted  here,  incorporates  a  digital  filter  to  predict  and  remove  the  acoustic 
signal.  We  design  this  digital  filter  using  the  same  system  ID  method  described  in  the  earlier 
section,  i.e.  using  the  system  ID  method  to  predict  the  acoustic  signal  with  the  actuators  on  and 
the  flow  off  and  then  subtracting  the  computed  acoustic  component  from  the  sensor 
measurement  with  both  the  actuators  and  flow  on  (see  Figure  5-6).  The  ID  parameters  were  p=l, 
n=100,  p  =  l.  Note  this  filter  has  much  higher  order  than  that  used  in  Figure  5-3.  However,  this 


74 


will  not  add  significant  computational  intensity  because  after  the  filer  is  designed  by  the  system 
ID  algorithm  the  filter  parameters  are  fixed  during  the  closed-loop  control.  Figure  5-7  shows  a 
comparison  between  the  actual  measured  acoustic  noise  with  wind  tunnel  off  and  the  predicted 
acoustic  noise  by  the  digital  filter.  Good  agreement  is  achieved.  To  test  this  digital  filter  further, 
we  used  the  same  digital  filter  in  Figure  5-7  but  reduced  the  actuator  amplitude  by  50%.  Figure 
5-8  shows  the  digital  filter  works  well  when  the  input  signal  is  changed,  indicating  the  linear 
behavior  of  the  acoustic  signal  produced  by  the  actuator  for  typical  excitation  levels  used  in  the 
experiment. 

Next,  the  digital  filter  is  applied  to  the  measurements  with  the  wind  tunnel  running. 
Figure  5-9  shows  the  power  spectra  comparison  of  the  pressure  measurements  before  and  after 
applying  the  digital  filter  for  a  system  ID  case.  The  power  spectrum  with  the  digital  filter 
applied  clearly  shows  lower  power  at  the  frequency  band  500  Hz  to  1500  Hz,  where  the 
piezoelectric  actuator  generates  most  of  its  acoustic  noise  due  to  the  actuator  resonance  as 
indicated  in  Figure  5-6.  The  digital  filter  is  thus  able  to  mitigate  the  acoustic  noise  component 
and  is  used  for  all  results  presented  below. 

5.1.3  Disturbance  Rejection 

5. 1.3.1  Closed-loop  control 

As  described  earlier,  the  disturbance  rejection  algorithm  requires  both  reference  (used  for 
feedback)  and  performance  signal  measurements.  As  described  earlier,  the  goal  of  the  algorithm 
is  to  minimize  the  fluctuations  of  the  performance  signal.  According  to  Venugopal  and 
Bernstein  (2000),  the  reference  and  performance  signals  can  be  the  same.  Herein,  different 
reference  and  performance  transducer  combinations  are  tested  for  comparison.  Since  SI  is  the 
closest  to  the  leading  edge  and  S6  is  the  closest  to  the  trailing  edge,  it  is  reasonable  to  investigate 
the  extremes  shown  in  Table  5-1. 

The  disturbance  rejection  algorithm  is  applied  to  all  four  cases  and  the  ID  and  controller 


parameters  are  summarized  in 
Case  Reference 

Performance 

cL 

cD 

L/D 

# 

Baseline 

transducer  y 

transducer  z 

0.21  ±0.02 

0.21  +0.09 

1.01  +0.08 

1 

SI 

SI 

0.84  +  0.01 

0.12  + 

6.97  +  0.37 

2 

S6 

S6 

0.83+0.01 

0.01 

0.12  +  0.01 

7.21  +0.46 

3 

SI 

S6 

0.84  +  0.01 

0.12  +  0.01 

7.11  +0.40 

4 

S6 

SI 

0.84  +  0.01 

0.12  +  0.01 

7.09  +  0.43 

Table  5-2.  Note  that  higher  values  increase  the  number  of  adjustable  parameters  in  the 
system  ID  and  the  controller,  which  may  increase  performance.  However,  this  is  not  guaranteed. 

First,  to  make  sure  the  control  is  inducing  a  global  effect,  we  examine  the  lift/drag 
performance  that  is  the  aerodynamic  objective  of  the  separation  control  system.  The  control 
objective  is  to  attach  the  separated  flow  and  thereby  reduce  the  fluctuating  pressure  spectrum 
associated  with  the  convection  of  the  vertical  structures  over  the  airfoil  surface.  The  lift  and  drag 
are  measured  after  the  closed-loop  control  algorithm  converges  by  the  balance,  which  is  only 
capable  of  providing  mean  or  time-averaged  data.  The  lift-to-drag  ratios  for  all  the  cases  are  also 
summarized  in  Table  5-1  and  include  uncertainty  estimates  that  account  for  calibration  and 
random  errors  (Tian  2007).  All  four  closed-loop  control  cases  give  similar  L/D  improvement, 
~  7x  L/D  of  the  uncontrolled  baseline  case.  Note  that  the  lift  is  increased  while  the  drag  is 
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decreased.  Close  inspection  via  tuft  and  smoke  flow  visualization  reveals  that,  in  all  4  cases,  the 
controller  is  able  to  fully  attach  the  separated  flow.  Since  the  four  closed-loop  control  cases  give 
the  same  L/D  within  experimental  uncertainty,  this  indicates  that  the  choice  of  the  performance 
and  reference  sensor  locations  does  not  have  a  significant  impact  on  the  integrated  lift/drag 
performance  for  this  flow  condition.  In  the  following  sections,  we  choose  to  study  case  #2  using 
S6  more  closely. 

5. 1.3.2  Effect  of  control  on  surface  pressure  signals 

Recall  that  we  use  the  actuator  and  surface  pressure  signals  to  model  the  plant  dynamics, 
and  the  disturbance  rejection  algorithm  attempts  to  suppress  the  surface  pressure  spectra.  Figure 
5-10  shows  the  time  traces  of  the  performance  surface  pressure  (measured  by  S6)  and  control 
input  signals  for  case  #2  before  and  after  the  control  is  turned  on.  The  results  clearly  show  that 
before  the  closed-loop  control  is  initiated,  the  performance  pressure  signal  has  relatively  large 
amplitude.  After  the  ID  and  control  is  initiated  from  a  zero  initial  condition  for  all  parameters, 
the  performance  pressure  signal  starts  to  decrease  driven  by  the  control  input  from  the  ZNMF 
actuators  generated  by  the  disturbance  rejection  algorithm.  The  entire  process  takes 
approximately  79  convective  time  scales  to  learn  the  dynamics  and  optimize  the  controller. 
After  a  steady  state  is  achieved,  the  pressure  signal  stays  at  the  lower  level  corresponding  to  an 
attached  low  oscillation  flow,  as  will  be  shown  in  the  next  section. 

Figure  5-1 1  shows  the  comparison  of  the  power  spectra  for  baseline  and  case  #2 
measured  by  the  performance  transducer  S6.  For  the  closed-loop  control  case,  the  power  spectra 
are  based  on  the  surface  pressure  signal  after  the  disturbance  rejection  algorithm  converges.  The 
noise  floor  of  the  transducer  is  also  plotted  for  comparison  to  verify  that  the  pressure  signals  for 
all  cases  are  well  above  the  noise  floor.  Figure  5-1 1  clearly  shows  that  the  disturbance  rejection 
algorithm  is  able  to  lower  the  spectrum  of  the  surface  pressure  signal  compared  with  the  baseline 
case  at  all  frequencies. 

5. 1.3.3  Quantitative  flow  visualization 

Preliminary  experimentation  shows  that  A1  and  A2  give  similar  results  but  A3  and  A4 
are  ineffective  because  they  are  located  downstream  of  the  separation  location.  Thus,  the  results 
with  A1  are  studied  closed.  Normalized  streamwise  velocity  and  vorticity  contours  (obtained 
using  500  PIV  image  pairs)  over  the  airfoil  for  the  baseline  and  closed-loop  control  case  #2  are 
shown  in  Figure  5-12  and  Figure  5-13.  For  the  closed-loop  control  case,  the  images  are  taken 
after  the  disturbance  rejection  algorithm  converges.  The  actuator  A 1  and  pressure  sensors  S 1  - 
S6  are  shown  as  circles  on  the  airfoil  surface.  For  the  baseline  case  in  Figure  5-12  (a),  the  flow 
separates  from  the  leading  edge  just  downstream  of  actuator  A1  and  all  six  pressure  sensors  are 
located  inside  the  separated  region.  Instantaneous  PIV  data  reveal  that  the  separated  flow 
features  large  coherent  vortices  sweeping  over  the  airfoil  upper  surface,  which  results  in  highly 
unsteady  pressure  signals  on  the  airfoil  upper  surface.  The  disturbance  rejection  algorithm 
“senses”  the  pressure  fluctuations  and  generates  actuation  signals  to  negate  the  pressure 
fluctuations.  This  process  ultimately  organizes  the  unsteady  flow  into  an  attached  turbulent  flow 
in  a  closed-loop  (smart)  fashion.  As  shown  in  Figure  5-12  (b),  the  flow  is  fully  attached  for  the 
closed-loop  control  case.  This  may  explain  why  the  four  cases  in  Table  5-1  result  in  similar 
lift/drag  performance,  i.e.  they  share  similar  information  about  the  flow  before  and  after  the 
closed-loop  control  is  initiated. 

5. 1.3.4  Control  input 

In  order  to  gain  physical  insight  into  the  actuator  control  signal  that  the  disturbance 
rejection  algorithm  generates  to  attach  the  flow,  we  examine  the  input  voltage  power  as  well  as 
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the  corresponding  input  electrical  power  to  the  actuator  Al.  The  first  quantity  has  units  of  V2, 
while  the  second  quantity  has  an  SI  unit  of  W.  Note  that  these  two  quantities  are  different 
because  piezoelectric  actuators  are  capacitive  devices  that  draw  higher  current  as  frequency  is 
increased  by  virtue  of  the  derivative  operation,  i=CdV/dt,  where  C  denotes  the  capacitance  of 
the  piezoelectric  actuators.  The  electrical  power  is  calculated  by  multiplying  the  input  voltage 
by  the  current  of  the  actuator.  The  current  is  measured  by  a  stand-alone  current  probe  (Tektronix 
TCP  A 300). 

The  voltage  and  rms  electrical  power  spectra  are  shown  in  Figure  5-14.  The  total  rms 
electrical  power  sums  up  to  12.7  mW.  It  is  clear  from  Figure  5-14  (a)  that  the  disturbance 
rejection  algorithm  generates  a  broadband  control  input  to  the  actuator  with  spectral  peaks  in 
both  the  low  (20  Hz  -  80  Hz)  and  the  high  (around  1  kHz)  frequency  ranges.  The  emphasis  at 
low  frequencies  is  due  to  the  larger  scale  coherent  flow  structures  described  earlier,  while  the 
emphasis  at  higher  frequencies  corresponds  to  the  smaller  scale  shear  layer  structures.  More 
detailed  discussions  on  the  two  types  of  flow  structures  can  be  found  in  Tian  et  al.  (2006)  and 
Wu  et  al.  (1998).  Since  the  flow  is  most  receptive  at  these  inherent  frequency  scales,  the 
disturbance  rejection  algorithm  attempts  to  utilize  the  two  characteristic  frequency  scales  by 
energy  addition  at  these  frequencies. 

On  the  other  hand,  actuator  characterization  experiments  reveal  that  the  actuator  produces 
very  small  output  at  low  frequencies  (<  500  Hz)  (Tian  et  al.  2006).  This  is  clear  from  Figure  5-5 
and  may  be  deduced  from  Figure  5-14  (b).  The  electrical  power  is  concentrated  at  higher 
frequencies  near  the  actuator  resonance  (than  the  voltage  signal  power)  while  remaining  almost 
flat  at  low  frequencies.  This  implies  that  although  the  controller  attempts  to  control  the  low  and 
high  frequency  instabilities  associated  with  the  wake  and  shear  layer,  respectively,  the  dynamic 
response  of  the  actuator  significantly  influences  the  control  system  dynamics.  It  is  important  to 
recall  that  the  “plant”  includes  the  actuator. 

Inspection  of  Figure  5-14  also  brings  into  question  whether  the  voltage  power  should  be 
used  as  a  penalty  function.  The  electrical  power  may  be  a  better  choice  for  the  penalty  function 
since  it  reflects  the  actual  power  consumption  by  the  actuator.  Such  questions  must  await  future 
studies. 

5. 1.3.5  Discussion 

To  examine  how  the  adaptive  controller  performs  under  different  flow  conditions,  the 
angle  of  attack  is  varied  continuously  from  12°  to  20°,  with  a  fixed  free  stream  Re  of  120,000. 
The  lift-to-drag  ratios  for  the  baseline  and  controlled  cases  are  plotted  in  Figure  5-15.  Clearly, 
the  baseline  flow  separates  starting  from  12°.  The  adaptive  controller  gives  the  best  performance 
at  AoA=12°.  The  performance  deteriorates  as  AoA  increases.  The  improvement  in  lift-to-drag 
ratio  becomes  very  small  at  AoA=20°,  which  means  the  controller  is  ineffective  at  this  AoA. 

The  comparison  of  the  power  spectra  for  baseline  and  case  #2  measured  by  the 
performance  transducer  S6  is  shown  in  Figure  5-16.  It  is  clear  that  unlike  the  12°  counterpart 
(shown  in  Figure  5-11),  the  power  spectrum  of  the  closed-loop  control  case  is  higher  than  that  of 
the  baseline  case.  This  defeats  the  purpose  of  the  closed-loop  controller  and  therefore  makes  the 
controller  ineffective.  The  key  question  is  why  the  pressure  spectral  characteristics  look 
conversely  different  at  different  AoAs.  Our  hypothesis  is  that  the  flow  can  only  be  partially 
attached  at  higher  AoA  in  the  mean  sense.  Instantaneously  the  partially  attached  flow  contains 
highly  unsteady  flow  vortices  that  cause  the  increase  in  the  pressure  spectra. 

Possible  solutions  include  using  alternate  surface  sensors,  such  as  MEMS-based  direct 
shear  stress  sensors  for  feedback  instead  of  pressure  sensors.  This  can  also  solve  the  acoustic 
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contamination  issue.  However,  we  believe  that  the  key  limitation  of  the  present  scheme,  as 
evidenced  by  its  failure  at  higher  angles  of  attack,  include  the  assumption  of  linearity  in  the 
system  identification  and  disturbance  rejection  algorithm.  The  nonlinearity,  especially  in  the 
partially  attached  unsteady  flow,  has  strong  impact  on  the  performance  of  the  closed-loop 
controller. 

Future  improvements  to  the  current  approach  include  exploration  of  nonlinear  control 
methods,  which  is  the  subject  of  the  next  section.  By  using  the  modulated  signals  as  input  (Tian 
et  al.  2006),  the  nonlinear  interactions  in  the  unsteady  flow  is  promoted.  On  the  other  hand, 
nonlinear  dynamical  approaches  are  desirable  compared  to  the  quasi-static  approach  in  our 
companion  study.  The  nonlinear  system  identification  algorithms  are  demonstrated  in 
simulations  by  Pillarisetti  and  Cattafesta  (2001).  Future  direction  includes  implementing  a 
nonlinear  dynamical  system  identification  and  control  scheme  in  the  wind  tunnel  experiments. 

5.2  Nonlinear  Control 

5.2.1  Experimental  Configuration 

As  shown  in  Figure  5-17,  dual-timing  control  loops  are  configured  to  implement  the 
optimization  algorithms  (described  in  the  next  section).  The  first  loop  synchronously  controls 
the  actuators  and  measures  the  low-pass  filtered  and  amplified  balance  signal,  while  the  second 
loop  averages  the  balance  output  and  performs  optimization  in  an  asynchronous  fashion.  The 
second  loop  acts  as  a  supervisory  controller  that  updates  the  control  parameters  in  the  first  loop. 
The  sampling  rate  of  the  first  loop  is  40  kHz,  while  the  second  loop  runs  on  a  host  PC  at  O(Hz). 
The  optimization  algorithm  is  programmed  in  Matlab  and  communicates  directly  with  the 
dSPACE  system  to  adjust  the  actuator  signal  parameters. 

5.2.2  Flow  Instabilities 

In  the  present  flow  conditions  (AoA=20°  and  Ref  =120,000),  the  baseline  uncontrolled 

flow  is  massively  separated  and  does  not  reattached  before  the  trailing  edge  (i.e.  post  stall).  As 
mentioned  earlier,  this  type  of  post-stall  flow  is  characterized  by  leading-edge  shear  layer  rollup 
and  vortex  shedding  in  the  wake  (Wu  et  al.  1998).  These  two  types  of  flow  structures  are  clearly 
visible  in  the  instantaneous  snapshot  of  the  flow  shown  in  Figure  5-18.  The  shear  layer  rollup 
structures  in  the  left  figure  have  a  much  smaller  length  scale  than  the  vortex  shedding  structures 
in  the  wake  shown  in  the  right  figure. 

To  characterize  the  velocity  fluctuations  in  the  two  flow  structures,  a  hot-wire 
anemometer  is  used.  The  hot-wire  is  traversed  vertically  across  the  shear  layer  (near  the 
separation  point)  and  the  wake  vortices  (1  chord  aft  of  the  trailing  edge).  The  maximal  u^ 

location  in  the  attached  sub-regions  is  then  determined  at  the  two  streamwise  locations.  Figure 
5-19  shows  a  plot  of  the  power  spectral  density  (PSD)  of  the  wake  and  the  shear  layer  at  the 
respective  peak  rms  locations.  The  PSD  was  estimated  using  a  4096  point  FFT,  a  Hanning 
window  with  75%  overlap,  and  320  effective  blocks.  The  plots  zoom  in  on  two  interesting 
regions.  The  left  plot  clearly  shows  the  dominant  wake  frequency  fwaJ((.  ~  40  Hz.  The  right  plot 

shows  the  much  higher  shear  layer  frequency  at  fSL~2040  Hz.  The  plot  also  provides  evidence 
for  the  nonlinear  coupling  between  the  shear  layer  and  wake  instabilities  via  the  presence  of  the 
shear  layer  frequency  fSL  and  the  sum/difference  frequencies  fSL±fwaXe-  N°te  that  fSL  is  much 
higher  than  fwake  in  accordance  with  classical  scaling  arguments  that  fsL~u~/escP  and 
fwake~u»/Wwakt  •  Based  on  the  definition  F+  =  fc/U„,  f=fSL  gives  F+~O(30)  and  f=fwake  give 
F+  ~O(0.6).  This  evidence  supports  our  hypothesis  that  more  than  a  single  characteristic 
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frequency  exists,  perhaps  explaining  the  wide  range  of  effective  forcing  frequencies  reported  in 
the  literature. 

To  study  our  hypothesis  about  the  nonlinear  quadratic  coupling  between  the  instabilities, 
higher-order  spectral  analysis  of  the  same  velocity  data  in  Figure  5-19  is  performed  (Nikias  and 
Petropulu  1993).  The  auto-bicoherence  contour  plot  shown  in  Figure  5-20  is  bound  between  0 
and  1  and  is  only  nonzero  due  to  nonlinear  quadratic  phase  coupling  (lock-on).  The  auto¬ 
bicoherence  thus  quantifies  the  fraction  of  power  in  a  random  signal  as  a  function  of  triad 
between  two  frequency  components  f, ,  f2  and  their  sum  or  difference  f,±f2 .  A  close  inspection 
of  the  contour  plot  reveals  distinct  features  at  f2=fw  (in  particular  near  fSL/2,fSL,1.5fsl ),  f2=fSL 
and  f 2  =fSL -f w  (in  particular  between  fSL  and  1.5fSL),  and  especially  along  the  lines  f,+f2=fSL  and 
f|+f2=fSL-fw .  These  data  conclusively  show  the  presence  of  nonlinear  quadratic  coupling 
between  the  Kelvin-Helmholtz  and  wake  instabilities. 

5.2.3  Actuator  Calibration 

5.2.3. 1  Frequency  response 

ZNMF  actuator  dynamics  is  a  critical  issue  for  the  control  of  a  separated  flow.  A  typical 
ZNMF  device  contains  a  cavity  and  a  vibrating  diaphragm  to  drive  oscillatory  flow  through  a 
small  orifice  on  the  cavity.  The  synthetic  jet  represents  a  coupled  electro-mechanical-acoustic 
system  with  frequency  dependent  properties  determined  by  device  dimensions  and  material 
properties  (Gallas  et  al.  2003).  Gallas  et  al.  use  the  lumped  element  modeling  approach  to  model 
these  actuators.  In  this  research,  although  modeling  our  ZNMF  devices  is  not  needed,  we  do 
need  to  characterize  the  frequency  response  of  the  actuator  Al  (that  is  used  in  this  research). 
Since  the  ZNMF  actuator  is  an  inherently  nonlinear  device,  the  traditional  approach  that  uses  a 
swept  sine  as  an  input  signal  is  not  appropriate.  Instead,  a  single  sine  wave  is  used  as  input,  and 
the  anemometer  signal  is  recorded.  The  frequency  is  then  increased  in  a  loop,  while  the  forcing 
amplitude  is  held  constant.  Figure  5-21  shows  the  rms  velocity  per  input  voltage  in  the 
frequency  band  from  500  Hz  to  2500  Hz.  Three  different  input  levels  are  used.  The  peak  output 
occurs  at  approximately  1200  Hz,  and  significant  output  is  apparently  limited  to  a  bandwidth  of 
500-1500  Hz.  The  output  level  is  very  low  for  frequencies  less  than  500  Hz  and  larger  than  1500 
Hz.  This  precludes  the  possibility  of  directly  forcing  either  the  low-frequency  (~40  Hz)  wake  or 
high  frequency  (-2020  Hz)  shear  layer  instabilities  via  sinusoidal  excitation.  It  also  highlights 
the  preferential  output  of  the  actuator  near  its  resonance  frequency.  Furthermore,  the  nonlinear 
nature  of  the  actuators  is  revealed,  since  the  frequency  response  function  is  not  independent  of 
the  input  voltage.  (If  the  actuator  were  linear,  these  curves  would  collapse.)  However,  that  is 
this  nonlinear  behavior  that  is  leveraged  to  enable  forcing  at  low  and  high  frequencies,  as 
explained  below. 

5.2.3.2  Types  of  actuation  waveforms 

Three  typical  multi-modal  waveforms  are  studied  to  take  advantage  of  the  multiple 
instabilities  of  a  separated  flow.  They  are  shown  in  Figure  5-22:  (a)  amplitude,  (b)  burst  and  (c) 
pulse  modulation.  In  (a)  and  (b),  the  lower  plot  in  the  figure  is  the  result  of  a  point-by-point 
product  of  the  top  two  waveforms.  Such  forcing  is,  in  general,  a  modulation  of  a  (usually)  high 
frequency  carrier  signal,  (e.g.,  a  sine  wave  with  frequency  fc )  by  a  low  frequency  modulation 
signal  (either  a  sine  wave  or  square  pulse  with  frequency  fm).  In  addition,  a  parameter  A  is 
multiplied  to  determine  the  amplitude.  For  the  BM  signal,  there  is  an  additional  parameter,  the 
duty  cycle,  which  determines  how  many  sine  wave  periods  occur  in  each  burst.  In  this  research, 
the  duty  cycle  is  adjusted  such  that  only  one  period  occurs  in  each  burst.  The  last  case  in  part  (c) 
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is  a  pulse  train,  which  can  be  interpreted  as  the  modulation  of  a  constant  signal  by  square  pulse. 
Similar  to  the  BM  signal,  the  duty  cycle  can  be  an  additional  parameter.  In  this  research,  it  is 
kept  to  be  the  shortest  achievable  width  on  the  dSPACE  control,  i.e.  one  discrete  sample  At .  In 
this  case,  there  is  only  one  waveform  parameter  fm  to  vary.  As  one  moves  from  amplitude  to 
burst  to  pulse  modulation,  the  modulation  process  results  in  an  increasingly  rich  signal  spectrum 
with  broader  spectral  content,  which  improves  the  likelihood  that  the  excitation  waveform  will 
excite  an  inherent  instability.  Furthermore,  the  required  actuator  power  reduces  as  will  be  shown 
in  the  following  section. 

Next  we  study  the  output  of  the  actuator  subject  to  the  AM  waveform  excitation  as  an 
example.  Figure  5-23  shows  the  velocity  output  of  the  ZNMF  actuator  A 1  measured  by  the  hot¬ 
wire  anemometry  subject  to  an  AM  excitation,  where  A=50  Vpp  (peak-to-peak  voltage),  fm=50 

Hz  and  fc = 1 1 80  Hz.  The  wire  was  placed  at  a  sufficient  distance  (~lmm)  above  the  actuator 

slot  so  that  there  is  no  reverse  flow  occurs  and,  hence,  no  signal  rectification  is  needed.  The 
temporal  record  shown  in  (a)  clearly  shows  the  low  frequency  oscillations  at  50  Hz  and  the  high 
frequency  oscillations  at  1180  Hz.  The  corresponding  power  spectral  density  of  the  velocity 
signal  shown  in  (b)  reveals  that  the  high  power  is  at  the  low  modulation  frequency  50  Hz  while 
the  second  highest  power  is  at  the  high  carrier  frequency,  with  additional  harmonic  distortion 
peaks.  Recall  that  the  response  to  the  sinusoidal  excitation  at  frequencies  less  than  500  Hz  is 
very  low,  as  shown  in  Figure  5-21.  The  modulation  enables  the  actuator  to  generate  high  power 
signals  at  low  frequencies,  while  the  sinusoidal  response  is  limited  by  the  actuator  dynamics. 
This  allows  an  actuator  operating  at  or  near  its  resonant  frequency  via  a  carrier  signal  at  fc  to 
generate  significant  disturbances  at  characteristic  frequencies  of  the  flow  that  are  far  from  the 
natural  frequency  of  the  device.  This  characteristic  is  attributed  to  the  nonlinear  nature  of  the 
actuator  system.  Similar  behavior  is  observed  for  BM  and  PM. 

5.2.3.3  and  electrical  power  calibration 

Figure  5-21  showed  that  the  actuator  ums  response  varies  significantly  with  actuation 

frequency  when  subjected  to  sinusoidal  excitation.  This  is  typical  for  a  ZNMF  actuator.  When 
the  performance  of  separation  control  is  studied,  the  actuation  frequency  is  always  of  prime 
importance.  However,  when  the  actuation  frequency  is  varied,  the  actuator  response  is  also 
varied  even  for  constant  amplitude.  As  mentioned  earlier,  it  has  been  shown  that  sinusoidal 
authority  varies  monotonically  with  Vj/U^  up  to  some  maximum  value  (Seifert  et  al.  1993, 
1996,  1999;  Glezer  and  Amitay  2002;  Mittal  and  Rampunggoon  2002).  In  other  words,  when  the 
actuation  frequency  is  varied,  the  two  performance-determining  parameters  (namely  frequency 
and  amplitude)  are  varied  simultaneously.  Unfortunately,  this  can  lead  to  misleading  results. 
For  example,  one  may  find  that  the  performance  is  the  best  at  the  peak  response  frequency  of  the 
actuator,  which  could  be  simply  because  the  actuator  is  providing  higher  output  as  opposed  to  the 
flow  being  more  receptive  as  that  frequency. 

The  same  problem  above  still  exists  even  when  multi-modal  waveforms  are  used;  the 
actuator  response  varies  when  A ,  fm  and  fc  are  varied.  To  separate  amplitude  forcing  effects, 
the  following  calibration  is  performed.  For  the  AM  and  BM  signals,  a  two-dimensional  grid  in 
the  (fm,fc)  space  is  generated,  and  CM  and  rms  electrical  power  consumed  by  the  actuator  are 

measured  for  various  excitation  amplitudes.  This  time-consuming  task  takes  days  to  complete 
and  verify  repeatability.  The  results  are  shown  in  Figure  5-24  and  Figure  5-25.  The  profiles  for 
five  amplitudes  (30  Vpp,  35  Vpp,  40  Vpp,  45  Vpp  and  50  Vpp)  are  recorded  and  shown.  The  lowest 
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profile  is  for  the  30  Vpp  case  and  the  highest  profile  is  for  the  50  Vpp  case.  For  the  PM  signal, 
there  is  only  one  frequency  parameter  fm;  thus  a  contour  plot  in  the  (A,fm)  space  is  shown  in 

Figure  5-26.  The  amplitudes  are  again  30  Vpp  to  50  Vpp  with  5  Vpp  increment. 

There  are  several  important  observations  in  the  results.  For  the  AM  signal,  the  response 
(CM  and  electrical  power)  is  approximately  independent  of  the  modulation  frequency  fm  but 

strongly  dependent  on  fc .  The  CM  at  a  fixed  fm  is  similar  to  the  sinusoidal  velocity  response 

shown  in  Figure  5-21.  In  addition,  the  shape  of  the  electrical  power  surfaces  are  somewhat 
different  than  the  shapes  of  the  ,  and  the  peak  frequencies  are  different  too.  This  means  that 

when  the  actuator  provides  the  maximum  velocity  output,  the  electrical  power  consumption  is 
not  maximized. 

For  the  BM  signal,  the  CM  and  electrical  power  responses  are  dependent  on  both  fm  and 
fc .  They  monotonically  increase  with  fm  as  expected.  As  with  the  AM  signal,  the  shapes  of  the 
CM  and  electrical  power  surfaces  are  different  and  possess  different  peak  frequencies.  The 
profiles  are  simpler  for  the  PM  signal  since  there  is  only  one  frequency  parameter.  Both  the  CM 
and  electrical  power  responses  increase  with  fm  and  A  as  expected.  Furthermore,  typical  levels 
of  CM  are  10~4 ,  1CT5  and  KT6  for  the  AM,  BM  and  PM  signals,  respectively. 

With  this  information  on  CM  and  electrical  power,  additional  constraint  functionality  is 
added  to  the  adaptive  optimization  program  to  hold  CM  or  electrical  power  constant,  while  fm 
and  fc  are  varied.  This  is  done  by  adjusting  A  in  accordance  with  Figure  5-24,  Figure  5-25  and 
Figure  5-26  in  each  iteration  of  the  optimization  algorithm.  Thus,  it  is  possible  that  at  certain 
(fm,fc)  combinations,  constant  CM  or  electrical  power  cannot  be  achieved.  In  this  case,  fm  and 

fc  are  set  to  the  values  at  the  boundary  such  that  the  momentum  or  power  constraint  is 

maintained.  This  will  hopefully  become  clearer  when  the  optimization  results  are  discussed  in 
the  next  section. 

5.2.4  Adaptive  Control  Results 

First,  the  constrained  optimization  is  carried  out  at  AoA=12“  and  Re=  120,000.  At  this 
AoA,  the  dynamic  control  approach  in  our  companion  study  works  well.  The  constraint  is  set  at 
CM  =  7.15  x  10  6.  This  is  a  relatively  low  level,  which  can  be  achieved  in  most  portions  of  the 

(fm,fc )  space  for  the  AM  and  BM  signals  and  all  values  of  fm  for  the  PM  signal.  To  protect  the 

actuator  from  physical  damage,  the  maximum  amplitude  is  limited  to  50  Vpp.  The  results  are 
summarized  in 


Table  5-3.  The  lift-to-drag  ratios  are  comparable  with  the  results  using  the  dynamic 
control  approach  (within  uncertainties).  This  is  not  surprising  because  the  flow  is  completely 
attached  by  both  approaches  at  AoA=12\ 

On  the  other  hand,  when  the  AoA  is  20“  where  the  dynamic  control  approach  fails  to 
attach  the  flow,  the  present  nonlinear  control  shows  its  clear  performance  improvement. 
Therefore,  a  more  detailed  study  is  performed  at  this  AoA. 
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In  addition  to  the  constrained  cases  for  =  7.15  x  10'6,  the  constrained  optimization  is 

also  carried  out  with  constant  electrical  power  of  0.0005  W.  The  constrained  optimization 
results  are  summarized  in  Table  5-4.  Typical  search  paths  are  plotted  in  Figure  5-27  to  Figure 
5-29.  The  shaded  areas  in  Figure  5-27  and  Figure  5-28  denote  the  achievable  parameter  space, 
inside  which  the  optimization  is  constrained.  It  is  clear  that  the  achievable  areas  cover  most  of 
the  parameter  space  for  the  AM  cases  while  limited  space  for  the  BM  cases.  Since  there  are  two 
parameters:  fm  and  fc,  three  initial  points  are  needed  in  each  optimization  run  as  mentioned  in 

earlier  sections.  The  symbols  A,  ■,  •  denote  the  search  paths  following  each  initial  point  during 
the  optimization.  On  the  other  hand,  for  the  PM  signal,  only  two  initial  points  are  needed. 
Different  sets  of  initial  points  are  chosen  to  cover  most  of  the  operational  space  to  achieve  a 
global  optimum. 

In  Table  5-4,  the  first  column  indicates  if  the  constraint  is  active  in  the  optimization  process. 

The  last  three  columns  summarize  the  optimal  values  for  converged  results.  The 
experimental  uncertainties  are  also  added  to  the  L/D  values.  By  simply  observing 
the  LTD  values,  we  find  that  the  performance  for  the  AM  case  (a)  gives  the  best 
performance.  The  AM  and  BM  signal  are  superior  to  the  PM  signal  with  constant  CM 

=  7.15  x  10  6.  For  the  AM  and  BM  cases  (a),  note  that  the  converged  fc  is  near  the 
shear  layer  frequency  fSL .  On  the  other  hand,  the  modulation  frequency  fm  assumes 
a  value  near  the  wake  frequency  fwake  and  its  superharmonics.  Most  importantly, 
effective  separation  control  is  achieved  by  using  the  multimodal  waveforms  with  CM 

of  at  least  an  order-of-magnitude  smaller  than  typical  values  reported  in  the  literature 
(Table  2  in  Greenblatt  and  Wygnanski  2000).  Selected  sinusoidal  excitations  are  also 
tested  on  our  airfoil  model  to  compare  with  the  results  using  the  multimodal 
waveforms.  The  voltage  amplitude  for  these  tests  is  held  at  50  Vpp.  The  results  are 
summarized  in 


Table  5-5.  It  is  clear  that  sinusoidal  control  (even  with  higher  CM)  gives  poorer 

performance  than  the  AM  and  BM  excitations. 

In  addition,  except  the  PM  signal,  the  nonlinear  control  is  able  to  achieve  much  better 
performance  than  the  dynamic  control  in  our  companion  study  at  AoA=20".  The  nonlinear 
control  approach  benefits  from  the  nonlinear  coupling  of  the  instabilities  and  the  integrated 
performance  (L/D)  measurements  instead  of  local  unsteady  pressure  measurements. 


Table  5-1.  Case  descriptions  and  performance  for  disturbance  rejection  experiments.  (AoA=12‘ 
and  Re=  120,000) 


Case  #  Reference 

transducer  y 

Baseline 

1  SI 


Performance 
transducer  z 

SI 


CL  CD  L/D 


0.21+  0.02  0.2 1  +  0.09  1.01+  0.08 

0.84  +  0.01  0.12  +0.01  6.97  +  0.37 
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2 

S6 

S6 

0.83  +  0.01 

0.12  +  0.01 

7.21  +0.46 

3 

SI 

S6 

0.84  +  0.01 

0.12  +  0.01 

7.11  +0.40 

4 

S6 

SI 

0.84  +  0.01 

0.12  +  0.01 

7.09  +  0.43 

Table  5-2.  Summary  of  parameters  in  disturbance  rejection  algorithm. 


P 

n 

P 

Pc 

nc 

Pc 

1 

2 

10 

p+n+p-1 

2 

20 

Table  5-3.  Constrained  optimization  results  using  the  AM,  BM  and  PM  signals.  (Baseline 
L/D=1.01  at  AoA=12"  and  Re=  120,000) 


Signal  type 

Constraint 

Constraint 

Active? 

Converged 

Converged 

f« 

Converged 

L/D 

AM 

C„  =  7.15  x  10'6 

No 

74 

1667 

7.47  +  0.45 

BM 

CM  =  7.15  x  106 

Yes 

48 

1305 

7.63  +  0.36 

PM 

C„  =  7.15  x  10'6 

No 

16 

NA 

7.14  +  0.25 

Table  5-4.  Constrained  optimization  results  using  the  AM,  BM  and  PM  signals.  (Baseline 
L/D=l.l  at  AoA=20°  and  Re=  120,000) 

Signal  type 

Constraint 

Constraint 

Active? 

Converged 

fm 

Converged 

fc 

Converged 

L/D 

AM:  Case  (a) 

C^  =  7.15  x  10'6 

Yes 

61 

2405 

2.18  +  0.07 

AM  :  Case  (b) 

Power=0.0005 

No 

202 

1005 

1.77  +  0.05 

BM  :  Case  (a) 

CM  =  7.15  x  10  6 

Yes 

55 

1979 

1.95  +  0.05 

BM  :  Case  (b) 

Power=0.0005 

Yes 

56 

1318 

1.52  +  0.04 

PM  :  Case  (a) 

CM  =7.15  x  10  6 

No 

16 

NA 

1.49  +  0.04 
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PM  :  Case  (b)  Power=0.005 


No 


29 


NA 


1.48  +  0.04 


Table  5-5.  Results  using  sinusoidal  excitations. 


F=fL.„/U, 

C, 

L/D 

0.5 

~0 

1.14  +  0.02 

15 

3.16  x  10  4 

1.76  +  0.03 

26 

7.79  x  10'6 

1.77  +  0.03 
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Dynamic  pressure 
transducers 
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Filter/ Amplifier 

Figure  5- 


Actuation 

signal 


(running  ID  and  control 
algorithms  in  real  time) 

I .  NACA  0025  airfoil  model  with  actuators,  sensors  and  closed-loop  control  system. 
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Figure  5-2.  Phase  averaged  pulse  response  measured  by  six  pressure  sensors.  The  slow 
propagation  velocity  of  the  coherent  flow  structures  is  clearly  visible. 
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Figure  5-3.  Comparison  between  measured  signal  from  the  pressure  sensor  (#6  in  Figure  5-1) 
and  the  fitted  output  by  ARMARKOV  system  ID  algorithm  for  long  and  short  time 
intervals.  Results  show  a  reasonable  match  at  low  frequencies  between  measured  and  fitted 
outputs.  For  ARMARKOV  ID:  p=l,  n=2.  Error!  Objects  cannot  be  created  from  editing 
field  codes.=  10. 


Figure  5-4.  Mean  Squared  Error  (Running  MSE)  between  measured  and  fitted  outputs.  Results 
show  that  the  ARMARKOV  ID  algorithm  converges,  i.e.  error  being  minimized. 
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Figure  5-5.  Comparison  between  frequency  response  (FR)  and  fitted  response  by  ARMARKOV 
ID  algorithm.  Parameters  for  FR:  Error!  Objects  cannot  be  created  from  editing  field 
codes.=4096  Hz,  NFFT=1024,  75%  overlap  and  Hanning  window.  For  ARMARKOV 
system  ID:  p=l,  n=2.  Error!  Objects  cannot  be  created  from  editing  field  codes.=10. 
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Figure  5-6.  Dual  signal  paths  from  the  actuator  to  the  pressure  sensor  (acoustic  and 
hydrodynamic).  A  digital  filter  is  introduced  to  remove  the  acoustic  component  by  turning 
off  the  flow  to  isolate  the  acoustic  path. 
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Figure  5-7.  Actual  measured  and  predicted  acoustic  Figure  5-8.  Actual  measured  and  predicted  acoustic 
noise  using  a  band-limited  random  signal  to  the  noise  using  the  same  filter  as  in  Figure  5-7  but 

actuator.  with  one  half  of  the  input  amplitude. 


Figure  5-9.  Power  spectra  of  the  sensor  signals  (with  wind  tunnel  running)  before  and  after 
applying  acoustic  filter. 
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Figure  5-10.  Performance  surface  pressure  (SI)  and 
control  input  signals  (in  Volt)  before  and  after 
the  ID  and  control  is  initiated  for  case  #2. 
Control  is  established  within  1  second  or  <100 
convective  time  scales. 


Figure  5-11.  Power  spectra  of  the  pressure  transducer 
output  for  the  baseline  and  the  closed-loop 
control  cases  measured  by  S6  (performance). 
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Figure  5-12.  Contours  of  streamwise  velocity  Error!  Objects  cannot  be  created  from  editing 
field  codes,  for  (a)  baseline  and  (b)  closed-loop  control  case  #2  at  AoA  =  12°  and 
Rec=  120,000. 
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Figure  5-13.  Contours  of  vorticity  for  (a)  baseline  and  (b)  closed-loop  control  case  #2  at  AoA  = 
12°  and  Rec=  120,000. 


Figure  5-14.  (a)  Voltage  and  (b)  electrical  power  spectra  of  the  actuator  A1  input  signal  for  the 
closed-loop  control  case. 
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Figure  5-15.  Performance  comparison  at  different  Figure  5-16.  Power  spectra  of  the  pressure  transducer 
AoA.  output  for  the  baseline  and  the  closed-loop 

control  cases  measured  by  S6  (performance)  at 
AoA  =  20°  and  Rec=  120,000. 
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Figure  5-18.  Flow  structures  in  separated  flow. 
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Figure  5-1 9.  Wake  (1  chord  aft  of  TE)  and  shear  layer  (near  separation) 
functions  at  peak  rms  location. 
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Figure  5-20.  Auto-bicoherence  of  the  same  velocity  signal  analyzed  in  Figure  5-19  using  the 
same  parameter  settings.  The  auto-bicoherence  is  zero  except  where  nonlinear  phase 
quadratic  phase  coupling  occurs  due  to  interactions  between  the  shear  layer  and  wake 
instabilities. 


Figure  5-2 1 .  Frequency  response  of  ZNMF  actuator  A 1 . 
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(a)  amplitude  modulation  (AM)  (b)  burst  modulation  (BM)  (c)  pulse  modulation  (PM) 

Figure  5-22.  Various  waveforms  of  unit  amplitude  A  =  1  that  can  be  used  to  excite  multiple 
instabilities  or  modes  in  a  separated  flow. 
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Figure  5-23.  Velocity  response  (a)  and  its  power  spectral  density  (b)  subject  to  an  AM  excitation 
for  the  ZNMF  actuator  Al.  A  =50  Vpp  ,  /m=50  Hz  and  /f  =  1180  Hz. 
Measurements  were  made  outside  the  region  of  reverse  flow. 
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Figure  5-24.  CM  (a)  and  electrical  power  (b)  profile  subject  to  AM  excitation.  Five  actuation 
amplitudes  (30  Vpp,  35  Vpp,  40  Vpp,  45  Vpp  and  50  Vpp)  are  shown. 
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Figure  5-25.  CM  (a)  and  electrical  power  (b)  profile  subject  to  BM  excitation, 
amplitudes  (30  Vpp,  35  Vpp,  40  Vpp,  45  Vpp  and  50  Vpp)  are  shown. 
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Figure  5-26.  CM  (a)  and  electrical  power  (b)  profile  subject  to  PM  excitation.  Five  actuation 
amplitudes  (30  Vpp,  35  Vpp,  40  Vpp,  45  Vpp  and  50  Vpp)  are  shown. 


Figure  5-27.  Constrained  search  using  AM:  Cases  (a)  and  (b).  The  shaded  area  denotes  the 
constraint  area  and  amplitude  is  adjusted  at  each  step  to  satisfy  the  constraint  criteria.  AOA 
=  20  deg.  and  Re  =  120,000. 


96 


Constrained  search 


Constrained  search 


Figure  5-28.  Constrained  search  using  BM:  Cases  (a)  and  (b).  The  shaded  area  denotes  the 
constraint  area  and  amplitude  is  adjusted  at  each  step  to  satisfy  the  constraint  criteria.  AOA 
=  20  deg.  and  Re  =  120,000. 
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6  Summary  and  Future  Work 


An  adaptive  system  identification  and  feedback  control  algorithm  is  applied  to  the 
separation  control  problem  for  a  NACA  0025  airfoil  at  nominal  angles  of  attack  of  12°  and  20° 
and  a  chord  Reynolds  number  of  120,000  with  a  tripped  boundary  layer,  corresponding  to  control 
of  a  massively  leading-edge  separated  flow.  In  particular,  a  recursive  ARMARKOV  system  ID 
algorithm  is  used  to  model  the  flow  dynamics  and  provide  the  information  required  to  implement 
the  disturbance  rejection  algorithm  in  real  time  with  no  prior  knowledge  of  the  system  dynamics. 
Phase-locked  PIV  and  fluctuating  surface  pressure  measurements  provided  evidence  of  the  link 
between  the  separated  flow  vertical  structures  and  the  surface  pressure  fluctuations.  The  chosen 
control  objective  was  thus  to  suppress  the  airfoil  surface  pressure  fluctuations.  The  disturbance 
rejection  algorithm  was  able  to  automatically  generate  control  input  to  the  ZNMF  actuator, 
emphasizing  low  (i.e.,  wake)  and  high  (i.e.,  shear  layer)  characteristic  frequencies  of  the 
separated  flow.  The  effect  of  the  control  is  to  enhance  near-wall  mixing  and  suppress  the  highly 
unsteady  flow  structures.  This  adaptive  control  scheme  is  able  to  completely  reattach  the  flow 
using  low  (-12.7  mW)  power  to  a  single  piezoelectric  synthetic  jet  actuator.  The  closed-loop 
control  results  show  -  7  x  improvements  in  the  lift/drag  ratio,  with  a  corresponding  increase  in 
lift  and  reduced  drag  and  concomitant  reductions  in  the  fluctuating  surface  pressure  spectra.  The 
present  results  are,  to  the  best  of  our  knowledge,  the  first  experimental  demonstration  of  an 
adaptive  dynamic  feedback  control  of  a  separated  flow.  The  results  reveal  the  tremendous 
potential  of  closed-loop  flow  control  to  real  aircraft  applications  but  also  reveal  key  issues 
worthy  of  further  study. 

First,  in  terms  of  positives,  the  adaptive  closed-loop  control  scheme  has  several  attractive 
features.  It  is  quite  general,  and  no  prior  knowledge  of  the  system  dynamics  is  required.  The 
system  identification  and  disturbance  rejection  algorithms  are  integrated,  and  the  system 
dynamics  can  be  obtained  with  minimal  a  priori  user  knowledge.  The  controller  is  implemented 
using  DSP  hardware  and  can  be  easily  incorporated  in  hardware-in-the-loop  applications.  It  can 
be  applied  to  not  only  flow  separation  control  problems,  but  also,  for  example,  cavity  oscillation 
control  and  turbulent  boundary  layer  control. 

Second,  in  terms  of  unresolved  technical  issues,  there  remain  many  concerning  the  actuator 
and  sensor  dynamics.  Clearly,  better  actuators  with  not  just  higher  output  but  flatter  dynamic 
response  over  wider  frequency  range  are  desirable.  While  measuring  unsteady  surface  pressure 
is  relatively  straightforward,  the  potential  acoustic  contamination  issue  was  highlighted.  This 
difficulty  was  mitigated  with  an  acoustic  digital  filter  here  but  at  the  cost  of  additional 
computational  complexity  from  an  already  limited  DSP.  This  issue  suggests  the  use  of  alternate 
surface  sensors,  such  as  MEMS-based  direct  shear  stress  sensors  for  feedback  instead  of  pressure 
sensors.  Or  perhaps  thermal  sensors  may  be  sufficient  despite  their  sensitivity  to  more  than  just 
shear  stress.  Ultimately,  it  is  believed  that  the  key  limitation  of  the  present  scheme,  as  evidenced 
by  its  failure  at  higher  angles  of  attack,  include  the  assumption  of  linearity  in  the  system 
identification  and  disturbance  rejection  algorithm.  The  second  part  explores  nonlinear  control 
methods. 

In  the  post-stall  separated  flow  where  the  flow  does  not  reattach,  there  are  two 
characteristic  instabilities:  the  shear  layer  Kelvin-Helmholtz  instability  and  wake  instability.  Our 
experiments  have  the  evidence  for  such  instabilities.  In  addition,  the  second  order  spectral 
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analysis  has  quantified  the  quadratic  phase  coupling  between  the  two  instabilities,  which 
indicates  the  separated  flow  is  a  complex  multi-frequency  system. 

Three  multi-modal  waveforms  (namely  amplitude  modulation,  burst  modulation  and 
pulse  modulation)  are  used  targeting  excitation  of  the  multi-frequency  separated  flow  system.  A 
simplex  optimization  approach  for  controlling  the  separated  flow  has  been  developed  to  search 
for  the  optimal  actuation  parameters  of  the  three  waveforms  using  the  ZNMF  devices.  It  is 
typical  for  the  CM  response  to  vary  when  the  waveform  frequency  parameters  vary  for  the 

ZNMF  devices.  This  can  potentially  lead  to  misleading  results  about  the  optimal  forcing 
frequency.  The  actuator  dynamics  is  taken  into  consideration  in  the  optimization  approach.  To 
offset  the  actuator  dynamics  implications,  a  special  routine  is  devised  to  hold  the  CM  at  constant 

during  the  optimization  process  utilizing  the  pre-calibrated  actuator  response  profiles.  This  is 
specifically  done  by  varying  the  actuation  voltages  according  to  the  response  profiles  to  keep  the 
CM  at  constant  levels. 

The  constrained  optimization  results  seeking  to  maximize  L/D  are  promising  and  reveal 
the  importance  of  forcing  nonlinear  interactions  between  the  shear  layer  and  wake  instabilities. 
Effective  separation  control  is  achieved  by  using  oscillatory  momentum  coefficients  o(l0^-105), 

which  is  more  than  an  order-of-magnitude  smaller  than  typical  values  reported  in  the  literature 
(see  summary  in  Greenblatt  and  Wygnanski  2000).  Specifically,  the  optimized  carrier  frequency 
fc  targets  the  shear  layer  frequency  while  the  optimized  modulation  frequency  fm  targets  the 
wake  frequency  and  its  super-harmonics.  The  nonlinear  control  is  able  to  achieve  similar 
performance  at  AoA=12”  and  much  better  performance  than  the  dynamic  control  in  our 
companion  study  at  AoA=20”.  The  nonlinear  control  approach  benefits  from  the  nonlinear 
coupling  of  the  flow  instabilities  and  the  integrated  performance  (L/D)  measurements  instead  of 
local  unsteady  pressure  measurements. 
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